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1.  INTRODUCTION 

The  advent  of  space-based  weapon  systems  in  our  times  has  raised  the  prospects  of  future 
'Star  Wars’  conflicts,  rendering  the  potential  use  of  explosive  devices  against  space  targets  a  present 
day  engineering  reality.  The  warhead  of  choice  in  space  seems  to  be  of  the  fragmentation  type,  for 
obvious  reasons.  The  effectiveness  of  fragments  is  unhampered  by  the  space  environment  (lack  of  air 
may  even  be  helpful).  By  contrast,  bare  charges  in  space  are  considerably  less  efficient  ttian  in  air. 
One  may  wonder  why  this  is  so  since  in  air,  as  in  space,  the  same  amount  of  chemical  energy  is 
released  through  the  detonation  process.  The  explanation  is  that  the  difference  is  in  the  much  larger 
mass  involved  in  the  air  blast,  relative  to  the  bare  charge  mass. 

For  a  more  comprehensive  explanation,  we  take  a  close  look  at  the  process  by  which  an  explosive- 
driven  air  blast  wave  is  generated.  The  explosive  products  effectively  constitute  a  rapidly  expanding 
spherical  piston  (typical  initial  speed  around  6  km/ sec),  which  drives  an  intense  shock  wave  into  the 
surrounding  air.  At  a  typical  range  of  100R0  (and  with  air  density  equal  to  about  1/1000  of  charge 
density),  the  mass  of  air  entrained  by  the  shock  is  about  1000  times  the  charge  mass.  Thus,  the 
highly  concentrated  initial  explosive  energy,  has  spread  over  a  much  larger  mass  than  that  of  the 
charge,  via  the  mechanism  of  wave  propagation  in  compressible  media,  resulting  in  an  increased 
momentum.  For  a  comprehensive  treatment  of  blast  waves  in  air  the  reader  is  referred  to  Baker[l]. 

It  is  also  worthwhile  noting  that  explosive  products  in  space  typically  attain  hypersonic  speed  prior 
to  impacting  at  the  target.  The  flow  velocity  in  an  air  blast  is  typically  subsonic  or  somewhat 
supersonic.  It  is  thus  expected  that  the  actual  gasdynamic  interaction  between  the  blast  flow  and  a 
stationary  target,  will  be  fundamentally  different  in  these  two  cases. 

We  contend  that  blast  effects  in  space  may  still  be  of  practical  interest  for  reasons  such  as  the 
following  : 

(i)  Notwithstanding  the  poor  efficiency  of  a  bare  charge,  its  use  should  not  be  ruled  out  altogether. 
Fragments  would  contribute  to  existing  -  and  potentially  hazardous  -  population  of  space 
debris,  underlining  the  obvious  fact  that  there  is  no  absolutely  safe  standoff  distance  from  an 
isotropic  fragmentation  warhead.  A  clean  bare  charge  may  thus  be  a  reasonable  alternative. 

(ii)  Even  a  fragmentation  warhead  has  some  residual  blast  capacity,  which  has  to  be  considered 
either  as  a  factor  in  enhancing  target  damage,  or  as  a  threat  to  be  reckoned  with  in  determining 
a  safe  standoff  distance. 
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The  key  idea  of  the  present  model  is  a  combination  of  the  assumption  that  target  dynamic 
response  is  related  primarily  to  total  blast  impulse,  and  the  physically  plausible  notion  that  this 
impulse  is  equal  to  the  total  momentum  of  that  portion  of  the  expanding  explosive  products  which 
impacts  at  the  target.  The  sense  in  which  this  simple  notion  constitutes  an  approximation  to  a  proper 
gasdynamic  analysis  of  the  interaction  between  the  fluid  and  the  target,  is  clarified  in  Ch.  2.  In  that 
chapter  we  also  present  an  illuminating  comparison  between  impulsive  blast  loading  in  air  and  in 
space. 

In  order  to  demonstrate  the  ChargeMass-Range-Damage  relationship  implied  by  our  impact  blast 
approximation,  we  chose  a  simple  target  model :  A  cantilever  beam  with  a  rigid-perfectly  plastic 
stress-strain  relationship.  It  represents  an  extended  structural  element  such  as  a  solar  panel  or  an 
antenna.  We  make  use  of  studies  conducted  by  Mentel  [21  and  by  Bodner  and  Symonds  [3J,  which 
showed  that  by  and  large,  the  effect  of  accelerating  the  beam  impulsively  was  to  cause  a  rotation 
about  a  plastic  hinge  at  the  point  of  support.  The  final  angle  of  rotation  is  generally  proportional  to 
the  initial  kinetic  energy,  so  that  equating  damage  with  that  angle,  results  in  damage  being 
proportional  to  the  square  of  the  impulse  imparted  to  the  target  by  the  blast  loading.  A  presentation 
of  this  dynamic  response  model,  including  a  sample  case,  is  given  in  Ch.  3. 

Our  ChargeMass-Range-Damage  relationship  may  imply  some  far-reaching  conclusions  when 
applied  to  the  analysis  of  a  more  general  configuration  than  the  single-charge/ single-target  case.  In 
Ch.  4  we  present  a  simple  analysis  of  a  sub-munition  configuration  of  N  bare  charges,  concluding 
that  it  seems  to  have  no  advantage  in  efficiency,  relative  to  a  single  charge  of  equal  mass.  Sections  5 
and  6  contain  conclusions  and  references,  correspondingly. 

We  conclude  the  introduction  by  listing  the  main  assumptions  made  in  the  present  study  : 

(a)  Blast  loading  and  target  response  are  uncoupled.  This  is  true  since  typically  the  target  mass  is 
much  larger  than  the  mass  of  that  portion  of  the  explosive  products  which  impacts  on  it. 

(b)  Dynamic  target  response  is  independent  of  specific  loading  time  history.  It  depends  solely  on 
total  (time-integrated)  impulse. 

(c)  The  target  is  a  panel  extended  as  a  relatively  supple  cantilever.  It  is  supported  by  a  relatively 
rigid  and  massive  core  structure. 

(d)  The  charge  is  a  sphere  detonated  at  its  center.  The  expansion  is  spherically  symmetric. 

(e)  Target  surface  is  normal  to  local  flow  vector. 

(f)  Target  orbital  velocity  relative  to  the  center  of  the  charge  is  negligible,  compared  with  the 
velocity  of  the  expanding  products. 
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2.  IMPACT  BLAST  LOADIN': 


Consider  the  expanding  explosive  products  impacting  at  a  target  as  shown  in  Fig.  2-1.  By 
regarding  the  fluid  as  an  ensemble  of  non-interacting  particles  moving  at  velocity  U(R,t) ,  and  by 
assuming  a  no- rebound  normal  impact  at  the  surface,  the  pressure  time  history  is  given  by  : 

P.(t)  -  p(R,<)|U(R,t)J2  (2-1) 


How  is  this  simple  impact  mechanism  related  to  the  actual  gasdynamic  interaction  between  the 
expanding  explosive  products  and  the  target?  When  a  target  is  located  at  a  range  of  at  least  several 
charge  radii,  two  features  in  the  free  stream  of  the  oncoming  fiuid  are  significant :  The  flow  is  highly 
hypersonic  (Mach  number  20  or  higher),  and  the  static  pressure  is  very  small,  which  means  that 
P  +  pU2  *  pU2  .  These  facts  were  bom  out  by  a  numerical  computation  which  we  performed  for  a 
typical  high  explosive  characterized  by  the  following  parameters : 

p0  ”  1800  (kg  m'3) 

Yd  ”  3 

(2-2) 

DCJ  -  8  (m  ms'1) 

Q0  -  dcj2/[2(Ycj2  "  01  -  4  (MJ  kg'1) 

Where  Q0  was  determined  by  assuming  that  the  detonation  corresponded  to  the  CJ  point  on  the 
explosive  Hugoniot  curve,  and  that  the  detonation  products  were  an  ideal  gas  with  a  specific-heat 
ratio  yCj  .  The  spherically  expanding  flow  was  computed  by  integrating  the  Euler  equations  for 
isentropic  flow  via  a  high-resolution  conservative  finite-difference  scheme  [4-6J.  The  initial  conditions 
were  the  self-similar  flow  field  of  a  just-detonated  spherical  charge  given  by  Taylor  [7].  The  code 
GRP  with  which  the  computation  was  performed  is  described  and  listed  in  Appendix  A. 

Consider  the  flow  at  a  stationary  target,  which  begins  at  the  moment  of  arrival  of  the  expanding 
explosive  products  (Fig.  2-2),  A  qualitative  description  of  the  ensuing  flow  pattern  is  made  by 
observing  its  evolution  in  time.  Immediately  following  the  initial  (normal)  impact,  the  fluid  is  stopped 
at  the  target  by  a  backward-propagating  shock  wave  reflected  from  the  surface.  Since  the  target  is  of 
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finite  extent,  the  fluid  between  the  shock  and  the  surface  is  accelerated  laterally,  and  streamlines  that 
tend  to  curve  around  the  target  are  being  formed.  If  the  oncoming  flow  were  stationary,  the  flow 
field  would  evolve  toward  the  familiar  configuration  of  a  detached  bow-shock  positioned  at  a 
relatively  narrow  standoff  distance  from  the  surface. 

Let  us  find  the  post-shock  pressure  in  these  two  limiting  phases.  In  the  initial  phase,  the  fluid  is 
stopped  at  the  target  by  a  reflected  shock  (Fig.  2-3a),  and  in  the  pseudo-stationary  phase  (Fig.  2-3b), 
the  shock  is  stationary.  In  either  case  we  find  the  post-shock  pressure  to  be  given  by  a  pressure- 
recovery  expression  of  the  form  : 

P2  -  apU2  (2-3) 

Where  a  is  a  consent  related  to  the  appropriate  y  (assuming  the  expanded  explosive  products 
are  an  ideal  gas).  The  governing  equations  in  the  reflected  shock  case  are  : 

p(U  +  S)  «  p2S 

p(U  +  S)2  ~  P2  (2-4) 

P(Y  +  1)/(Y  “  1)  “  p2  (strong  shock) 

Where  the  unknowns  are  p2  ,  P,  ,  S  . 

The  equations  for  the  stationary  shock  case  are  : 
pU  *  p2U2 

pU2  -  P2  +  P2U22  (2-5) 

P(Y  +  1)/(Y  ~  1)  =  p2  (strong  shock) 

Where  the  unknowns  are  p7  ,  U2  ,  P,  .  Thus,  solving  for  a  in  the  two  cases  represented  by 
equations  (2-4)  and  (2-5),  we  get  : 
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Reflected  shock  a.  ”  fty  +  l)/2]2 

(2-6) 

Stationary  shock  a  «*  2/(y  +  1) 

In  either  case,  since  the  gas  is  not  dense,  the  effective  range  cf  y  is  somewhere  between  1.0  and 
1.4  ,  so  that  setting  a  =  1  is  an  approximation  commensurate  with  the  overall  crudeness  of'  the 
present  impact  blast  model.  Since  the  flow  in  the  layer  between  the  shock  and  the  target  is  low 
subsonic  (at  least  it  is  so  away  from  target  edges),  the  post-shock  pressure  is  a  reasonable  substitute 
for  the  surface  pressure.  Also,  a  »  1  is  an  appropriate  approximation  where  the  flow  is  so 
rarefied  that  it  is  collisionless.  In  this  limit,  a  *  1  corresponds  to  full  thermal  accommodation  of 
re-emitted  molecules  from  a  presumably  cold  surface. 

The  foregoing  analysis  constitutes  a  justification  of  the  impact  approximation  to  the  surface 
pressure  (2-1).  Now  we  turn  to  the  task  of  evaluating  the  impulse  which  is  defined  as  the  time- 
integrated  surface  pressure.  Using  the  impact  approximation  (2-1),  the  impulse  is  given  by  : 

to  CO 

p,(t)dt  =  Jp(R,l)[U(R,t)l2dt  (2-7) 

o  o 

Let  us  introduce  a  Lagrange  mass  coordinate  m  which  enables  a  transformation  from  the  Euler 
system  (R,t)  to  the  Lagrange  system  (m,t).  The  differential  relation  associated  with  this 
transformation  at  constant  R  is  : 

dm  -  4;cR2p(R,t)U(R,t)dt  (2-8) 


Since  it  is  assumed  that  the  fluid  is  not  accelerated  at  any  (R,t)  in  the  range  of  interest  for  blast 
loading,  the  velocity  U(R,t)  can  be  regarded  as  function  solely  of  the  mass  coordinate  ,  so 
that  U(R,t)  =  U(m)  .  Using  (2-8)  we  are  then  able  to  cast  the  impact  blast  expression  (2-7)  in  the 
following  simple  and  physically  appealing  form  : 


I(R)  -  Z/4jtR2 

W 

Z  =  /  U(m)dm 


/■ 


(2-9) 


The  upper  limit  VV  in  (2-9),  which  is  consistent  with  the  upper  limit  °o  in  (2-7),  implies  that  the 
total  impulse  is  somewhat  overestimated,  since  it  contains  contributions  from  the  innermost  layers  of 
the  explosive  products  that  will  arrive  at  the  target  as  t  -»  oo  , 


The  total  momentum  Z  is  thus  a  constant  which  can  be  evaluated  for  any  specific  explosive 
charge  by  numerical  integration.  We  performed  this  computation  with  the  code  GRP  described  in 
Appendix  A.  In  doing  so  for  the  typical  explosive  (2-2),  we  found  out  that  the  impulse  (2-9)  was  a 
reasonable  approximation  at  ranges  as  low  as  R  -  3R0  .  Furthermore,  it  was  found  that  Z  could 
be  approximated  by  the  maximum  attainable  momentum  for  the  given  charge  mass  and 
energy  W(2Q<p1^2  ,  to  within  about  6%  .  Apparently,  the  total  momentum  is  not  overly  sensitive 
to  the  exact  velocity  distribution  function  U(m)  ,  so  that  assuming  a  value  of  Z  appropriate  to  the 
uniform  distribution  U(m)  •  (2Q^)'/2  is  a  reasonable  approximation.  Thus  we  finally  arrive  at  the 
following  closed-form  approximation  for  the  blast  impulse  : 


I(R)  -  K  W(2Q0),/1/4rtR2 

K  -  1 


(MO) 


Where  the  coefficient  K  is  retained  in  order  to  suggest  that  its  value  be  determined  more  accurately 
from  detailed  experimental  or  computational  data,  in  the  event  that  such  data  become  available.  At 
present  our  best  estimate  is  K  -  1 . 


There  is  one  companson.  however,  which  can  readily  be  made  with  available  data.  We  refer  to 
impulsive  blast  loading  in  air.  such  as  given  by  Baker  (Ref.  1,  Fig.  6.3  in  the  supplement).  The 
companson  is  conveniently  made  with  a  non-dimensional  form  of  (2-10),  which  is  rewritten  as  : 

I  -  !<R)[4kR02/W<2Q0),'2|  -  (R/R0)2  (2-11) 

The  air  blast  data  has  to  be  converted  to  the  same  normalizattoi  cheme  as  .  ■  E4.  (2-11),  before 
the  companson  can  be  made.  Considering  the  definition  of  1  in  (2-11)  above,  and  the  definition  of 
scaled  range  arid  air  blast  impulse  (Table  6.2  of  Ref.  i ),  this  conversion  is  done  by  multiplying  the 
sealed  an  impulse  and  range  by  the  following  coefficients  (sea  level  aii  is  issu  t  1) : 


Impulse  Multiplier 
Range  Multiplier 
P,  -  1.3  (kg  m  J) 


P  *  3(2y)*,/2(4rt/3)1/3  (P,/P,Q0)^6  (P,/P0)^2  -  .01204 
6  *  (4n/3)'/3  <P0Q0/P.)!/3  -  67.06 


(2-12) 


P#  *  0.1  (MPa) 
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Y  *  1.4 


The  air  blast  conversion  was  done  by  a  small  code  which  is  given  in  Appendix  B.  The  air  and 
space  blast  impulses  are  shown  in  Fig.  2-4.  We  note  that  at  ranges  larger  than  about  10  charge  radii, 
the  air  blast  impulse  is  higher  than  the  space  impulse,  and  the  gap  widens  as  the  range  increases. 
This  observation  is  consistent  with  the  qualitative  explanation  given  in  the  introduction,  which 
attributed  this  effect  to  the  increase  in  the  entrained  air  mass  at  higher  range.  At  ranges  lower  than 
10  charge  radii,  the  air  mass  is  relatively  insignificant,  so  that  one  may  expect  the  blast  impulses  in  air 
and  in  space  to  be  comparable  Indeed,  the  inverse-square  variation  of  impulse  with  range  is 
apparent  for  the  air  blast  at  low  range.  In  absolute  values,  however,  the  low-range  space  impulse  is 
higher  by  a  factor  of  about  1.7.  This  might  be  interpreted  as  indicating  that  choosing  k  «■  1/1.7 
would  be  the  appropriate  "calibration’.  However,  we  do  not  propose  to  do  so,  since  we  are  not  able 
to  trace  the  various  factors  affecting  the  low-range  impulse  as  given  by  Baker  [1];  they  may  somehow 
depend  on  the  presence  of  air,  as  well  as  on  other  parameters  such  as  target  size  and  equation  of  state 
of  the  explosion  products. 


Figure  2-1.  Impact  Blast  Loading 


REFLECTED  SHOCK 


Figure  2*2.  Shock  Reflection  at  Impact  Phase 


TARGET  TARGET 


(a)  Initially  Reflected  Shock  (Impact )  (b  )  Stationary  Shock 


Figure  2-3.  Limiting  Cases  of  Shock  Reflection 


Figure  2-4.  Impulse  of  Normally  Reflected  Blast  Wave  at  Sea- Level  and  in  Space 


3.  TARGET  DYNAMIC  RESPONSE 


For  the  sake  of  constructing  representative  ChargeMass-Range-Damage  relations  from  our  impact 
approximation  to  the  blast  impulse  (2*10),  we  suggest  a  simple  idealized  structure  as  target  model.  It 
is  a  cantilever  beam  made  of  a  metal  characterized  by  a  rigid-perfectly  plastic  stress-strain  relation. 


This  model  is  supposed  to  represent  an  extended  spacecraft  component  such  as  a  solar  panel  or  an 
antenna.  The  core  structure  is  assumed  to  be  much  more  massive  and  rigid  than  the  extended 
structural  element,  so  that  the  cantilever  can  be  idealized  as  being  rigidly  supported.  The  sole 
dynamic  and  structural  parameters  are  hence  those  of  the  cantilever. 

For  this  purpose  we  make  use  of  an  experimental  and  theoretical  investigation  of  uniform 
cantilever  beams  subjected  to  impulsive  loading  that  was  conducted  by  Mentel  [2].  Aluminum  alloy 
beams  were  held  in  a  massive  support  that  was  gliding  along  a  rail  at  speed  V  ,  until  it  was  abruptly 
stopped  by  a  very  massive  anvil.  After  the  system  came  to  rest;  the  beams  were  observed  to  have 
rotated  through  an  angle  ©  about  the  point  of  support,  with  little  deformation  elsewhere  (Fig.  3*1). 


fj 
1  | 

u 


The  theoretical  model  suggested  by  Mentel  [2]  for  predicting  ©(V) ,  can  be  described  as 
comprising  two  stages.  Immediately  following  the  impact,  the  beam  comm-nces  rotating  rigidly 
about  the  support  point,  with  an  angular  momentum  equal  to  the  pre-collision  moment  of 
momentum  about  that  point.  This  application  of  the  principle  of  conservation  of  moment  of 
momentum  entails  an  abrupt  re-distribution  of  velocity  in  the  beam,  with  velocity  being  proportional 
to  distance  from  support,  and  the  tip  moving  at  1.5  V  .  The  angle  0  is  subsequently  determined 
from  the  requirement  that  the  rotational  kinetic  energy  be  dissipated  as  plastic  hinge  work  Mp0  . 
The  resulting  0(V)  expression  is : 

0  -  (3/8)pLV2/Mp  (3-1) 


r 


We  now  make  one  more  step  in  formulating  the  model,  in  that  we  postulate  that  the  angle  0  is 
a  measure  of  damage  .  Using  the  following  expressions  for  Mp  ,  p  and  V  : 


Mp  -  (l/4)YIr 


M  =  Pph 


(3-2) 


V  -  I (R)/p 
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We  get  from  (2-10)  and  (3*1)  the  following  ChargeMass-Range-Damage  (W-R-0)  relationship  : 


R  «  CW1^ 

C  -  K3/167t2e)(LQ0/ppYh3)l1/4 


(3-3) 


We  note  that  the  effective  nnge  for  a  specified  target  and  "damage  level"  0  ,  is  proportional  to 
the  square  root  of  the  charge  mass  W  . 


Using  the  data  for  the  typical  explosive  (2-2),  and  the  following  data  for  a  specific  aluminum  beam, 
we  get  for  this  sample  case  : 


h  -  0.002  (m) 


L  -  1.0  (m) 

pp  -  2700  (kg  m (3-4) 
Y  -  300  (MPa) 

C  »  1.85  O'1'4  (m  kg'1/2) 


The  ChargeMass-Range-Damage  relationship  corresponding  to  this  sample  case  is  depicted  in 
Fig.  3-2. 
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4.  CLUSTER  CONFIGURATION 


In  a  cluster  configuration,  the  gain  in  damage  is  presumably  a  result  of  a  favorable  design  tradeoff 
between  reduced  charge  mass  and  reduced  range.  Can  such  a  gain  be  achieved  for  a  space  system, 
assuming  the  ChargeMass-Range-Damage  law  (3-3)  to  hold?  It  can  be  shown  that  by  adopting 
some  simple  strategy  of  sub-munition  dispersion  and  initiation,  equation  (3-3)  implies  no  gain  in 
target  damage. 

Let  us  assume  for  the  sake  of  a  reasonably  simple  analysis,  that  dispersion  and  initiation  of  sub¬ 
charges  would  take  place  according  to  the  following  scheme  : 


(a)  The  N  sub-charges  appear  to  fan  out  from  a  common  virtual  center,  moving  at  equal  speeds. 
At  subsequent  times,  their  centers  are  uniformly  distributed  over  an  expanding  spherical 
envelop. 

(b)  The  target  moves  at  a  constant  velocity  relative  to  the  virtual  center.  Its  point  of  closest 
approach  to  that  center  is  at  range  R  . 

(c)  The  timing  for  dispersion  is  chosen  so  that  the  target  intersects  (tangentially)  with  the  spherical 
envelop  at  the  point  of  closest  approach  (Fig.  4-1).  This  is  also  the  point  at  which  the  blast 
from  a  single-charge  configuration  detonated  at  the  virtual  center,  would  have  impacted  at  the 
target. 

(d)  All  sub-munitions  are  detonated  at  this  "moment  of  closest  approach". 

(e)  It  is  assumed  that  each  spherical  cap  of  area  4tiR2/N  will  contain  one,  and  only  one,  sub¬ 
charge.  The  probability  of  the  charge  location  on  that  cap  is  assumed  to  be  uniformly 
distributed.  The  expected  location  on  the  cap  is  hence  that  latitude  line  <p  which  divides  the 
cap  into  two  parts  of  equal  area  (Fig,  4-2). 

(f)  It  is  assumed  that  the  target  is  subjected  to  the  blast  of  a  single  sub-charge,  which  is  located  on 
the  mid-area  latitude  <p  of  the  spherical  cap  that  surrounds  the  target  (Fig.  4-2). 

Since  the  area  of  the  spherical  cap  subtended  by  <p  is  4nR2/(2N) ,  the  angle  <p  is  given  by  : 

sin(q>/2)  =  (2N)',/f2  (4-1) 

We  seek  a  comparison  between  the  deflection  0  for  a  single  charge  (W,R)  ,  and  the  deflection 

0^  in  the  sub-munition  case  (  =  W/N  v  R^  =  2Rsin(<p/2)  ).  From  the  ChargeMass-Range- 

Damage  law  (3-3),  using  also  Eq.  (4-1),  we  get : 
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(V«)  -  OVN/W)J  (R/RN)4  -  1/4 


(4-2) 


Consequently,  there  is  no  potential  gain  in  a  tradeoff  between  charge  mass  and  range,  for  a  cluster 
configuration  with  the  aforementioned  dispersion  scheme.  The  factor  1/4  ,  along  with  the  mass 
overhead  inherent  in  constructing  a  multi-charge  configuration,  indicate  that  in  causing  blast  damage, 
a  single  charge  is  more  effective  than  an  equal-mass  isotropically  dispersed  cluster. 
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5.  DISCUSSION  AND  CONCLUSIONS 


Our  analysis  pertains  to  a  bare  explosive  charge  initiated  at  a  point  of  closest  approach  to  the 
target.  We  have  shown  that  the  loading  impulse  on  a  planf'orm  target  is  given  by  the  impact 
approximation  (2*7),  which  states  that  the  impulse  is  proportional  to  the  charge  mass  and  inversely 
proportional  to  the  range  squared.  The  impulse  in  space  has  been  compared  with  impulse  in  air  at 
sea-level.  It  was  found  that  the  two  are  quite  comparable  at  close  range  (10  charge  radii  or  less), 
exhibiting  identical  variation  with  range.  At  far  ranges,  the  impulse  in  air  is  the  higher  one.  This  is 
consistent  with  the  notion  that  spreading  the  explosive  energy  over  larger  air  mass  results  in  larger 
momentum  (and  hence  reflected  impulse).  We  then  proceeded  to  develop  the  ChargeM ass- Range- 
Damage  law  (3-3)  for  an  impulse-responsive  target,  w’hich  states  that  blast  damage  is  proportional  to 
the  square  of  the  charge  mass  and  inversely  proportional  to  the  fourth  power  of  the  range.  These 
results  were  obtained  by  introducing  extensive  simplifications  in  the  analysis  of  gasdynamic 
interaction,  and  in  the  analysis  of  dynamic  target  response.  We  have  further  shown  that  this  damage 
law  also  implies  that  no  gain  can  be  achieved  by  an  idealized  cluster  configuration  of  bare  sub¬ 
charges,  relative  to  a  single  charge  of  equal  total  mass. 

It  is  worthwhile  noting  that  all  assumptions  introduced  in  the  course  of  formulating  the  impact 
blast  approximation  and  the  structural  dynamic  response  to  impulsive  loading,  imply  that  target 
damage  is  overestimated.  The  only  exception  is  the  approximation  in  setting  «  =»  1  ,  which  can  be 
readily  rectified  by  assigning  to  a  the  reflected  shock  value  given  in  (2-6).  Furthermore,  we  assumed 
that  the  pressure  at  the  midpoint  of  the  target,  is  the  pressure  everywhere  on  the  target.  Due  to  flow 
around  the  edges,  the  average  pressure  is  lower  than  the  midpoint  pressure.  Also,  targets  are  not 
everywhere  normal  to  the  flow  (and  charge/ target  attitude  is  not  a  design  parameter).  Oblique  impact 
obviously  entails  reduced  target  loading.  In  the  area  of  structural  dynamic  response,  a  time- 
distributed  loading  function  generally  delivers  less  kinetic  energy  to  the  structure  than  an  impulsive 
loading  of  equal  total  impulse,  resulting  in  reduced  deformation  (damage).  Thus,  while  the  present 
model  may  be  regarded  as  an  over  estimate  when  applied  to  a  sure-fail  analysis,  it  is  particularly 
suitable  in  determining  a  sure-safe  range. 
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APPENDIX  A.  The  GRP  Code 


The  purpose  of  this  Appendix  is  to  provide  a  concise  description  of  the  GRP  code,  and  a  listing  of 
its  CHARGE  version.  It  is  intended  for  users  that  have  had  prior  experience  in  implementing 
schemes  for  solving  the  Euler  equation  of  compressible  flow.  The  theoretical  background  of  GRP 
schemes  constitutes  the  principles  on  which  the  code  is  founded.  Some  familiarity  (at  least)  with  this 
background,  as  given  in  References  4,  S  and  6  is  indispensable  to  any  implementation  of  GRP 
schemes.  Reference  4  is  recommended  as  an  introduction.  The  planar  GRP  scheme  is  fully  described 
in  Reference  5,  and  the  duct-flow  GRP  scheme  on  which  the  present  CHARGE  version  is  based  is 
given  in  Reference  6.  (In  CHARGE  version  the  flow  is  spherical  and  the  “duct"  area  is  set  to 
X(I)**2  ,  but  the  code  can  handle  any  area  variation  -  see  subroutines  CROSS  and  RATIO  below). 

In  GRP  schemes,  second-order  accuracy  is  achieved  by  considering  a  piecewise  linear  interpolation 
of  the  flow  in  each  cell  (Fig.  A-l),  from  which  second-order  accurate  fluxes  at  each  cell  interface  are 
evaluated  through  an  analysis  of  a  local  Generalized  Riemann  Problem  (GRP).  Briefly  stated,  the 
GRP  goes  one  step  further  than  the  Riemann  Problem  (RP),  in  that  it  seeks  (analytically)  the  first 
time-derivative  of  the  flow  that  evolves  as  the  "diaphragm"  is  removed  from  the  cell  interface,  at  the 
origin  of  the  centered  (X,T)  wave  paths  of  the  RP  solution.  The  major  computational  subroutines  arc 
CYCEUL  where  the  integration  of  conservation  laws  is  performed,  RIEMAN  where  the  local 
Riemann  Problems  are  solved  by  Newton- Raphson  iterations,  MAGA  where  the  closed-form 
expressions  derived  from  the  GRP  analysis  [6]  are  used  to  compute  flow  time-derivatives  along  the 
contact  surface,  FLUXE  where  ail  the  previously  computed  information  is  used  to  extrapolate  the 
fluxes  to  mid-time-step  (T+DT/2)  which  constitutes  a  second-order  accurate  flux. 

The  plan  of  tnis  Appendix  is  as  follows.  Array  variables,  including  those  which  carry  conserved 
variables  (mass,  momentum  and  energy),  are  described  in  section  A.l.  This  is  followed  by 
descriptions  of  general  parameters  (A. 2),  labeled  COMMON  variables  (A. 3)  and  all  subroutines  (A. 4). 
We  conclude  by  giving  the  CHARGE  version  listing  (A. 5),  which  should  be  consulted  whenever  a 
reading  of  this  code  description  is  attempted. 

NOTE  :  The  present  CHARGE  version  was  implemented  in  a  GRP  code  version  that  had  been 
converted  to  treat  detonation,  weaves  as  chemically  reactive  compressible  flow.  However,  the 
detonation  scheme  is  effectively  neutralized  by  setting  QDET  =  0  (in  NETUNM).  All  variables 
pertaining  to  detonation,  such  as  arrays  Z(I),  DZ(I),  FIMZ(I),  ZMDOT(I)  and  labeled  COMMON 
variables  containing  Z  in  their  names,  should  be  ignored. 


A.!  Array  Variables 


The  code  GRP  is  organized  so  that  all  major  subroutines  are  called  with  standard  list  of  array 
variables  which  represent  the  integration  scheme  (i.e.  the  conservation  laws),  local  Riemann  Problem 
solutions  and  second-order  accurate  fluxes.  Virtually  all  array  variables  are  initially  defined  in 
BEGIN  (initial  conditions),  and  are  subsequently  updated  at  each  time  step  in  CYCEUL.  The 
following  list  explains  the  meaning  of  these  variables.  Some  terms  used  in  the  list  are  defined  below. 


X(I) 

U(l) 

P|I) 

RO(I) 

E(I) 


DU(I) 

DP(I) 

DRO(I) 

DG(I) 

DXSI(I) 

MIN(I) 

L’S(I) 


PS(I) 


L’IDOT(I) 

PIDOT(I) 

FIMZ(I) 

ZMDOT(I) 


grid  point  coordinate, 
velocity  in  cell !. 

pressure  in  cell  I  (computed  from  equation  of  state). 

density  in  cell  i.  This  variable  is  time-mtegrated  according  to  the  law  of 
conservation  of  mass.  (Com!  uted  in  CYCEUL). 

total  energy  per  unit  volume  (including  kinetic  energy)  in  cell  I.  This  variable  is 

time -integrated  according  to  the  law  of  conservation  of  (total)  energy.  (Computed 

in  CYCEUL). 

velocity  difference  in  cell  L 

pressure  difference  in  cell  I. 

density  difference  in  cell  I. 

Lagrange  sound  velocity  difference  in  cell  I. 

the  Lagrange  coordinate  increment  defined  as  RO(I)*(X(I  +  l)-X(I)),  for  cell  I. 
inactive  in  present  version. 

velocity  at  the  contact  surface  obtained  after  the  resolution  of  the  local 

discontinuity  at  X(I)  (Riemann  Problem  solution).  It  is  denoted  as  L'1  in 

References  '4-6. 

pressure  at  the  contact  surface  obtained  after  the  resolution  of  the  local 

jJi 

discontinuity  at  X(I)  (Riemann  Problem  solution).  It  is  denoted  as  P  in 

References  4-6. 

time  derivative  of  US(I)  along  the  contact  surface.  (This  derivative  is  the  result  of 

the  GRP  analysis.  It  is  computed  in  MAGA.  See  Ref.  5  and  6). 

time  derivative  of  PS(I)  along  the  contact  surface.  (This  derivative  is  the  result  of 

the  GRP  analysis.  It  is  computed  in  MAGA.  See  Ref.  5  and  6). 

inactive  in  present  version. 

inactive  in  present  version. 
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TENA(I) 

FIRCKI) 

FIM(I) 

FIE(I) 

G1P(I) 

VOL(I) 

Z(I) 

DZ(I) 


momentum  per  unit  volume  RO(I)*U(I)  in  cell  I.  This  variable  is  time-integrated 
according  to  the  law  of  conservation  of  momentum.  (Computed  in  CYCEUL). 
mass  flux  at  point  X(I)  (second-order  accurate), 
momentum  flux  at  point  X(I)  (second-order  accurate), 
energy  flux  at  point  X(I)  (second-order  accurate). 

the  pressure  term  in  the  momentum  flux.  It  corresponds  to  G(U)  in  References  4 
and  6. 

volume  of  cell  I. 
inactive  in  present  version, 
inactive  in  present  version. 


Glossary  of  terms  used  in  the  array  variables  list : 

Cell  I  -  the  cell  between  grid  points  X(I)  and  X(I  + 1).  All  cell  variables  are  averages  per  that 
interval. 

Difference  in  cell  I  -  the  difference  between  values  of  variable  at  cell  boundaries  X(I  + 1)  and  X(I). 
Those  values  are  obtained  from  "monotonized"  piecewise  linear  distribution  of  each  variable 
in  each  cell.  (Fig.  A-l). 

Second-order  accurate  flux  -  the  flux  time-derivative  at  point  X(I)  is  computed  from  the  time- 
derivatives  of  pressure  and  velocity  along  contact  surface  PIDOT(I)  and  UIDOT(I)  (in 
FLUXE).  Then  the  the  flux  is  extrapolated  to  the  centered  time  point  (T  +  DT/2),  using 
those  derivatives.  This  centered  value  is  the  second-order  flux  for  integrating  the 
conservation  laws  between  T  and  T  +  DT. 
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AJ  Major  Parameters 


A  list  of  major  parameters  indicating  their  meaning  end  the  routine  in  which  they  are  defined,  is 
given  below.  Those  parameters  defined  in  NETUNM  are  the  run  input,  There  is  no  reading  of  an 
input  file  in  this  version  of  GRP  code  (and  the  only  output  is  the  printed  output). 


L 

LL 

T 

DT 

TMAX 

TMUD 

DTMUD 

NCYC 

COLELA 

KEYMON 

NCYCPR 

STAB 

DTBA 

DTKOD 

KDT 


number  of  grid  points  +•  1  (main  program) 

L  -  1  (MAIN  PROGRAM) 
time  (MAINO) 
time  step  (MAINO) 

maximum  time  (when  T.GE.TMAX  the  run  is  terminated)  (NETUNM) 
time  for  which  next  printing  will  take  place  (NETUNM) 
printing  time  step  (NETUNM) 

serial  number  of  time  step  (integration  cycles)  (MAINO) 

switch  to  evaluate  cell  differences  by  Colella's  method  when  COLELA.NE.O 
(NETUNM) 

key  for  monotoruzation  scheme  (just  one  is  presently  provided  when  COLELA. EQ.O) 
(NETUNM) 

frequency  ofline  priming  at  each  cycle  (time  step)  (NETUNM) 

CFL  stability  coefficient.  Must  be  smaller  than  1.  (NETUNM) 
next  time  step  computed  from  stability  criterion  (CYCEUL) 
former  time  step  (MAINO) 

index  of  cell  where  DTBA  was  determined  (CYCEUL) 


A.3  Labeled  COMMON  reliables 


Labeled  COMMONS  are  used  primarily  to  transmit  data  to  and  from  routines  that  perform  the 
major  computational  steps  of  the  GRP  scheme,  .i.e,  RIEMAN,  MAGA  and  FLUXE;  these  routines 
are  called  from  CYCEUL,  When  the  value  of  any  of  those  variables  is  needed  for  later  use,  whether 
for  updating  conservation  variables  (RO,  TENA,  E),  or  for  printing,  it  is  stored  in  the  appropriate 
array.  All  labeled  COMMON  variables  are  grouped  under  labels  that  indicate  their  role,  and  their 
names  are  also  mnemonic.  Generally,  suffix  L  means  Left  and  suffix  R  means  Right.  It  may  indicate 
sides  either  with  respect  to  a  cell  interface  X(l),  or  with  respect  to  the  contact  surface  which  separates 
the  Right-  and  Left-  propagating  waves  in  a  solution  to  the  local  Riemann  Problem.  We  indicate  by 
INPUT  variables  tha*  are  computed  prior  to  calling  the  subroutine,  and  by  OUTPUT  variables  whose 
value  was  computed  within  the  subroutine  and  constitutes  the  result  of  calling  that,  subroutine. 

COMMON  /5TEP0/  Parameters  related  to  the  local  Riemann  Problem.  This  is  the  first  step  in 
the  GRP  scheme. 

UL,  PL,  ROL,  CL,  GL,  SL  -  velocity,  pressure,  density,  sound  speed,  Lagrange  sound  speed  and 
entropy,  attributed  to  Left  side  of  cell  interface  at  point  X(I).  (INPUT) 

U'STAR,  PSTAR  -  velocity  and.  pressure  at  the  contact  surface  obtained  when  the  local 
discontinuity  is  resolved  (i.e.,  the  solution  to  the  local  Riemann  Problem).  The 
omission  of  L  or  R  suffix  indicates  that  P  and  U  are  continuous  across  the  contact 
surface.  (OUTPUT) 

RSTARL,  CSTARL,  GSTARL  -  density,  sound  speed  and  Lagrange  sound  speed  on  the  Left  side  of 
the  contact  surface.  (OUTPUT) 

WL  -  Lagrange  velocity  of  propagation  of  the  Left-moving  shock,  relative  to  the  fluid.  (OUTPUT) 
UW(6)  -  velocity  of  propagation  of  each  wave  front  (Fig.  A-3),  relative  to  the  inertial  system 
(X).  (OUTPUT) 

HELEML  -  logical  variable.  If  HELEML.EQ..TRUE.  the  Left-propagating  wave  is  a  shock. 

Otherwise  it  is  a  (centered)  rarefaction  wave.  (OUTPUT) 

NFLUX  -  integer  variable.  It  denotes  the  region  in  the  Riemann  solution  wave  structure,  which 
contains  the  point  X(I)  lor  all  time.  Refer  to  Fig.  A-3  for  illustration.  (OUTPUT) 
LAMDAL,  RATEL,  TEMPL,  TEMPSL,  ZL,  ZSTARL  -  inactive. 


COMMON  /STEP!/  Parameters  related  to  the  time-derivative  evaluation  of  the  GRP  scheme, 
performed  in  MAGA.  The  time-derivatives  of  P  and  U  along  the  contact  surface  are 
the  main  result  of  MAGA. 

DU1DT,  DPIDT  -  time-derivatives  of  velocity  and  pressure  along  contact  surface.  (OUTPUT) 

ASTARL  -  The  directional  derivative  of  U  along  the  fan  characteristic  at  the  trailing  characteristic 
of  the  Left  rarefaction  wave.  It  is  not  evaluated  when  the  Left  wave  is  a  shock.  (See 
References  4-6)  (OUTPUT) 

DGIDTL,  DRIDTL  -  time-derivatives  of  Lagrange  sound  speed  and  density  along  the  left  side  of 
the  contact  surface.  (OUTPUT) 

DSDAL  -  Lagrange  spatial  derivative  of  entropy  on  the  left  side  of  contact  surface,  prior  to 
removal  of  the  partition  at  X(I). 

SH,  RAT  -  the  cross-section  area  and  the  x-derivative  of  ln(SH).  They  are  user-defined  in  CROSS 
and  RATIO  respectively. 

DSDASL  -  entropy  derivative  used  in  the  special  'sonic"  case  (i.e,  when  NFLUX-2  or 
NFLUX-  5).  See  References  5,6  for  details.  (OUTPUT) 

LAMDSL,  DZDAL,  BETACL,  DZDASL  -  inactive. 


COMMON  /GRADS/  Used  to  transmit  flow  gradients  (that  exist  in  fluid  prior  to  removal  of  the 
partition  at  X(I))  to  MAGA. 

DUDXIL,  DPDXIL,  DGDXIL,  DRDXIL,  DSDXIL  -  gradients  of  U,  P,  G,  RO,  S  (with  respect  to 
Lagrange  coordinate).  They  are  computed  in  CYCEUL  for  transmission  to  MAGA. 
(INPUT) 

DZDXIL  -  inactive. 


COMMON  /FI/  Used  to  return  values  of  updated  flux  and  cell-interface  variables  from 

FLUXE. 


FJH1,  FIH2,  FIH3  -  second-order  flux  of  mass,  momentum  flow  (just  RO*U'**2)  and  energy.  They 
are  extrapolated  to  Half  the  time  step  T  +-  DT/2.  (OUTPUT) 

GIH  -  the  value  of  P  at  T+  DT/2 


UXN,  PXN,  GXN,  ROXN  -  values  of  L\  V,  G,  RO  extrapolated  to  New  time  T+DT,  at  celt- 
interface.  They  are  used  in  CYCEUL  to  get  tentative  (prc-monotonized)  new  cell 
differences.  (OUTPUT) 

ZXN,  FIH4,  ZMDOTL,  ZMDOTR  •  inactive. 


A.4  Description  of  Subroutines 


MAIN  PROGRAM 

The  task  of  this  program  is  to  allocate  array  space  for  the  NMAT  arrays  required  by  the  present 
version  of  GRP  code.  The  length  of  each  array  is  L.  The  allocation  is  done  by  calling  MAINO.  This 
standard  calling  sequence  is  maintained  hereafter,  thus  facilitating  modifications. 

MAINO 

This  subroutine  functions  as  an  overall  organization  routine.  It  can  be  read  as  a  kind  of  flow¬ 
chart  of  the  entire  computation.  First,  run  set-up  is  done  by  calling  once  to  NETUNM  (data)  and 
BEGIN  (initial  conditions).  Then  a  loop  over  time  steps  is  begun.  In  each  cycle  the  integration  by 
one  time  step  is  performed  by  calling  CYCEUL,  and  subsequently  boundary  conditions  arc 
implemented  by  calling  SAFAE.  Whenever  T.EQ.TMUD.  results  are  printed  by  calling  PRINT  and 
TMUD  is  updated  by  adding  DTMUD. 

NETUNM 

Here  data  are  set  for  a  particular  run.  User  is  invited  to  modify  this  routine.  There  is  no  input 
file.  This  routine  is  called  just  once  from  MAINO.  Note  that  the  detonation  data  section  is  skipped 
when  QDET.EQ.O. 

BEGIN 

Initial  conditions  are  set-up  in  this  routine.  The  configuration  of  some  nominal  case  is  given  in 
present  version.  (In  CHARGE  version  it  is  the  detonated  spherical  charge,  using  the  Taylor  self 
similar  solution  as  initial  conditions).  User  is  called  to  modify  this  routme  so  as  to  generate  any  other 
desired  initial  configuration. 

TAYLOR 

The  purpose  of  this  routine,  along  with  ancillary  routines  INIDAT,  RUNGE  and  DER1V,  is  to 
compute  the  self-similar  Taylor  solution  [7]  of  a  detonated  spherical  charge,  and  implement  it  as 
initial  conditions  for  the  GRP  computation  of  the  ensuing  expansion.  TAYLOR  is  called  once  by 
BEGIN. 

The  core  of  the  solution  is  the  numerical  (Runge-Kutta)  integration  of  two  coupled  ordinary 
differential  equations.  The  integration  variable  is  PS!.  (The  flow  velocity  normalized  by  DCJ  is  given 
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by  U-EXP(-PSI) ).  The  two  dependent  variables  are  X  -  the  normalized  radial  coordinate  (X-  1  at 
the  sphere  boundary),  and  C  -  the  normalized  speed  of  sound.  The  integration  is  carried  out  by 
calling  RUNGE,  which  in  turn  calls  DER1V  for  the  evaluation  of  derivatives.  Data  for  the  TAYLOR 
computation  is  set  up  by  calling  (just  once)  INIDAT. 

The  initial  conditions  needed  in  BEGIN  are  values  of  mass,  momentum  and  (total)  energy  per  cell. 
These  are  most  accurately  computed  by  spatially  integrating  the  Taylor  solution,  resulting  in  lumped 
mass,  momentum  and  energy  per  cell,  which  are  then  divided  by  the  cell  volume.  This  refinement  is 
significant,  since  gradients  are  high  near  the  charge  boundary  (X-  1).  A  total  mass  and  energy  check 
for  the  entire  sphere  is  performed  and  printed. 


INIDAT,  RUNGE,  DERIV 

Subroutines  used  only  in  conjunction  with  the  Taylor  initial  conditions  setup.  See  TAYLOR 
above. 

RATIO,  CROSS 

User-defined  routines.  If  A(X)  is  the  duct  cross-section  area,  then  CROSS(X)-*  A(X)  and 
RATIO(X)=D[ln(A(X))]/DX  . 


CYCEUL 

This  is  the  central  computation  routine.  All  major  stages  of  the  GRP  scheme  are  performed  by 
calling  specific  subroutines  from  CYCEUL.  Then  RO(I),  TENA(I)  and  E(I)  are  updated  to  new  time 
T  +  DT  by  solving  the  appropriate  conservation  laws  in  CYCEUL. 


The  first  loop  (DO  1)  performs  a  set  of  preparatory  steps  as  follows  : 

(a)  CALL  RIEMAN  -  Solving  the  local  Riemann  Problem  at  each  X( I). 

(b)  CALL  MAGA  -  Solving  the  local  Generalized  Riemann  Problem  at  each  X(l). 

(c)  CALL  FLUXE  -  Computing  second-order  fluxes  at  Xy I). 

(d)  Evaluation  of  cell* interfact  finite  differences  DU(I),  DP(I),  DRO(I)  in  each  cell.  These  will  be 
used  at  the  future  time  step  (after  monotonization)  for  piecewise-linear  interpolation  of  the  flow 
in  each  cell.  (See  definition  of  DUDXIL,  DPDX1L,...,  just  preceding  the  call  to  MAGA  in  this 
loop). 


Note  that  in  present  CHARGE  version  additional  computation  of  PRESS,  PULSE1,...,  PULSE4 
has  been  added.  It  is  just  informative  and  does  not  interfere  in  any  way  with  the  execution  of  the 
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GRP  scheme.  The  purpose  of  this  computation  is  to  monitor  the  numerical  solution  and  to  observe 
the  accuracy  within  which  the  asymptotic  value  of  the  momentum  integral  Z  (Eq.  2-9  above)  is 
approached. 

In  the  second  loop  (DO  2),  the  integration  of  the  three  conservation  laws  is  performed,  using 
second-order  fluxes  that  had  been  computed  in  loop  1.  Flow  variables  such  as  P(I)  and  U(I)  are 
computed  in  this  loop  from  the  conserved  variables.  The  cycle  computation  is  concluded  by  calling 
BDOK1  for  monotonization  of  DU(I),  DP(I)  and  DRO(I). 

SAFAE 

In  this  routine  user-defined  boundary  conditions  are  implemented.  Present  version  (CHARGE) 
contains  rigid  wall  at  the  center  of  the  sphere  X(2)-0,  and  an  "open  boundary"  at  the  outer 
computational  zone  limit  X(L).  The  rigid  wall  condition  is  achieved  by  setting  up  a  virtual 
antisymmetric  cell  next  to  the  boundary  cell,  so  that  the  solution  to  the  local  Riemann  Problem  will 
result,  in  a  non-moving  contact  surface  (USTAR=*  0).  The  open  boundary  is  an  approximation  to  an 
ideally  non-reflecting  boundary.  Here  the  virtual  cell  is  I*L,  and  the  flow  in  it  is  defined  as  a 
"continuation"  of  the  flow  in  the  adjacent  last  ceil  I  =*  LL. 

BDOK1 

Here  the  tentative  cell- interface  differences  DV(I)  are  monotonized  according  to  neighboring 
average  cell  values  V(I-l),  V(I)  and  V(I+l).  The  basic  idea  is  that  the  cell-interface  slope  DV(I) 
should  have  the  same  sign  as  the  average  slope  V(I  +  l)-V(I-l).  When  V(I)  is  a  local  extremum  DV(I) 
is  set  to  zero.  Also,  the  absolute  value  of  DV(I)  is  constrained  so  that  the  jump  from  a  cell-interface 
value  to  the  adjacent  average  value  V(I),  will  never  be  of  opposite  sign  to  DV(I). 

DCOLE 

When  COLELA  option  is  used  (not  in  present  CHARGE  version),  the  pre-monotonized  slopes 
are  simply  the  centered  difference  (V(I  +  l)-V(I-l))/2.  Note  that  even  under  this  option,  the 
monotonization  routine  BDOK  l  is  subsequently  called. 

PRINT 

Printing  of  results.  Reading  this  routine  is  self-explanatory.  Note  some  features  added  for 
present  CHARGE  version.  User  is  called  to  modify  this  routine  to  his  specific  needs. 
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SOF 


Run  termination  when  an  error  has  been  detected.  I  STOP  is  an  informative  index.  Ail  printing 
of  relevant  information  should  be  done  at  the  calling  routine  prior  to  calling  SOF.  Note  that  the  run 
is  ended  in  SOF  by  deliberately  causing  a  system  error  of  computing  SQRT(-l).  This  is  done  in  order 
to  trigger  printing  of  the  sequence  of  calling  routines  by  the  operating  system. 


RIEMAN 


Here  a  single  Riemann  Problem  (RP)  is  solved  by  calling  RIEMAN  from  CYCEUL.  Referring 
to  Fig.  A-2,  the  RP  is  solved  by  finding  the  point  of  intersection  (USTAR.PSTAR)  of  Left- 
propagating  and  Right-propagating  shock/rarefaction  adiabats  in  the  (U,P)  plane.  Prior  to  the  actual 
computation,  the  qualitative  wave  structure  is  determined.  It  is  characterized  by  the  index  NCASE  as 
follows : 


NCASE*  1 
NCASE -2 
NCASE -3 
NCASE  *4 


Left  wave  is  rarefaction,  Right  wave  is  shock. 
Both  waves  are  shock. 

Left  wave  is  shock,  Right  wave  is  rarefaction. 
Both  waves  are  rarefaction. 


The  computation  of  (USTAR.PSTAR)  is  coded  separately  for  each  case.  Newton-Raphson 
iteration  is  employed,  the  first  guess  being  the  intersection  of  the  Left  and  Right  rarefaction  branches 
(or  their  extrapolations),  which  is  done  in  closed-form.  Since  in  a  smooth  flow  this  guess  is  close  to 
the  exact  (USTAR.PSTAR),  little  extra  CPU  effort  is  spent  on  subsequent  Newton-Raphson 
iterations.  These  are  truly  needed  only  in  regions  of  shock  wave  computation. 

The  computation  in  RIEMAN  is  concluded  by  computing  UW(1) . UW(5)  (UW(6)- infinity). 

From  these  wave  speeds,  the  flux  index  NFLUX  that  denotes  the  location  of  the  X-axis  on  the  (X,T) 
wave  diagram  of  the  RP  solution  (Fig.  A-3),  is  evaluated.  It  is  later  needed  in  subroutine  FLL'XE. 
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MAGA 


The  major  purpose  of  this  routine  is  to  compute  DUIDT  and  DPIDT  along  the  contact  surface 
of  the  RP  solution.  Since  U  and  P  are  continuous  across  the  contact,  so  are  their  time-derivatives 
along  the  contact.  Thus,  DUIDT  and  DPIDT  are  solved  from  a  set  of  two  linear  equations.  The 
coefficients  of  each  equation  are  determined  by  GRP  analysis  of  the  wave  on  one  side.  See 
References  4-6  (particularly  Ref.  6)  for  details. 

FLUXE 

The  major  task  of  this  routine  is  to  compute  second-order  fluxes.  This  is  done  in  two  phases. 
The  first  phase  is  up  to  statement  9  CONTINUE,  where  using  NFLUX  the  X-sufRxed  values  of  flow 
variables  and  their  time-derivatives  are  defined.  An  X-suffix  means  that  the  variable  or  its  time- 
derivative  are  related  to  the  line  X“X(I)  on  the  (X,T)  wave  diagram  (Fig.  A-3).  In  the  second  phase, 
these  variables  and  their  time-derivatives  are  used  to  extend  fluxes  at  X(I)  to  Half-time-step  (hence 
the  suffix  H),  i.e.  T  +  DT/2.  It  is  these  fluxes  which  are  the  second-order  accurate  fluxes  for  the 
integration  of  the  conservation  laws  from  T  to  T+DT.  Also,  cell-interface  flow  variables  (suffix  N) 
are  extended  to  New  time  level  T  +  DT.  These  are  later  used  in  defining  cell  differences  DU(I).  DP(I) 
and  DRO(I)  in  CYCEUL. 


A.5  Listing  of  GKP  Code 


v&Rstoti 


C60PTI0NS  LIST 

IMPLICIT  REALX8CA-H,0-Z,*> 

C  PROGRAM  GRP  -  GENERALIZED  RIEMANN  PROBLEM. 

C  EXPANSION  OF  A  DETONATED  SPHERICAL  CHARGE  IN  VACUUM. 

C  INITIAL  CONDITIONS  FROM  TAYLOR'S  SELF  SIMILAR  SOLUTION. 

COMMON  BC102,26) 

I  , ENDB 

COMMON  /AB/A<50) 

EQUIVALENCE  ( L , Ail ) ) , ( LL , A( 2) ) , CT, A(3) ) , ( DT, A< A) ) , < TMAX, ACS) ) , 

1  (TMUD,AC6)),CDTMUD,AC7)),CJOB,A(8>),CNERI,AC9)), 

2  CJJJ,AC1Q)), (KEYMON, AC  11) ) , CNCYC, AC  12) ) 

CCQLELA, AC  13) ) 

CLAGEUL, AC 14) ) 

( UGAL, AC  15) ) 

CKEY£K,AC16)) 

( NCYCPR, AC  17 ) ) 

CSTAB,AC18)),CDTBA,AC19)),  CDTK0D,AC20>), CKDT,AC21)) 
/M0NIT/CASEAVC4),NC14C4),NF16C6), 

NMONU( 4 ) , NMONP C  4 ) , NMONRO (4), NMONZ  <  4 ) 


EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
COMMON 


CHACOOl 

CHA0002 

CHA0003 

CHA0Q04 

CHA0005 

CHA0006 

CHA0007 

CHA0008 

CHA0009 

CHAOOin 

CHAOOll 

CHA0012 

CHA0013 

CHA0014 

CHA0C15 

CHAQ016 

CHA0017 

CHA0018 


CHA0019 

DIMENSION  NZEROC  26 )  CHA0020 

EQUIVALENCE  C NZEROC 1 ), NC14C 1 ) )  CHA0021 

COMMON/PULS/PRESSC 1 0) ,  PUI.SE1  CIO),  PULSE2C 10) ,  PUI.SE3C  10)  , PULSE4C 10)  CHA0022 
CXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA0023 


DO  20  N»l,26 

20  NZEROCN) =0 
DO  21  N=l,4 

21  CASEAVCN)=0. 

DO  31  N=l,10 
PRESSCN) =0 . 

PULSE1(N)=0 . 

PULSE2(N) =0 . 

PULSE3(N) =0 . 

PUL5E4CN)=0. 

31  CONTINUE 

NMAT=26 

)  L=CL0CFCENDB)-L0CFCBC1,1)))/NMAT 

L  =  1 02 
LL=L-1 
NN=NMATXL 
DO  1  J  =  1,L 
DO  1  11=1, NMAT 
1  B( I , II  )  =  0  . 

CALL  MAINOC  L ,  BC  1 »  1),BC1,  2),BC1,  3),BC1,  4),B(1,  5), 

1  BC 1 ,  6  ) , DC  1 ,  7 ) , BC 1 ,  8 ) , BC 1 ,  9),BC1,10), 

2  BC1,11),BC1,12),BC1,13),BC1,14),BC1,15), 

3  BC1,16),BC1,17),BC1,18),BC1,19),BC1,20), 

4  BC1,21),BC1,22),BC1,23),BC1,24),BC1,25), 

5  BC 1 , 26  )  ) 

STOP 

_ EMIL 


CHA0024 

CHA0025 

CHA0026 

CHA0027 

CHA0028 

CHA0029 

CHA0030 

CKA0031 

CHA0032 

CHA0033 

CHA0034 

CHA0035 

CHA0036 

CHA0037 

CHA0038 

CHA0039 

CHA0040 

CHA0041 

CHA0042 

CHA0043 

CHA0044 

CHA0045 

CHA0046 

CHA0047 

CHA0048 

CHA0049 

JIHAQQJiCL 


SUBROUTINE  MAINO  M/H/VO  CHA0051 

1  CL,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN,  CHA0052 

2  US, PS, UI DOT , PI DOT ,  CHA0053 

X  FIMZ, ZMDOT ,  CHA0054 

3  TENA, FIRO, FIM, FIE,GIP,VOL,Z, DZ)  CHA0055 

IMPLICIT  REALX8CA-H,Q-Z,$)  CHA0056 

DIMENSION  XCL),UCL),PCL),ROCL),GCL),ECL),DUCL),DP(L), DROC L ) ,  CHA0057 

1  DGCL),DXSICL),MINCL),  CHA0058 

2  USCL),PSCL),UIDOTCL),PIDOTCL)  CHA0059 

3  ,TENACL),FIROCL),FIMCt.),FIECL)  CHA0060 

4  ,GIPCL),VOLCL),ZCL),DZCL)  CHA0061 

5  , FIMZC  L ) , ZMDOT  CL)  CHA0062 

COMMON  /AB/AC50)  CHA0063 

EQUIVALENCE  C L L , AC  2 ) ) , C T , AC  3 ) ) , C DT , AC  4 ) ) , C TMAX , AC  5 ) ) ,  CHA0064 

1  C  TMUD, A(6 ) ) , C  DTMUD, AC  7 ) ) , C  JOB, AC8) ) , (NERI , AC  9 ) ) ,  CHA0065 

2  CJJJ.AC10) ),C  KEYMON, A  C 1 1 ) ) , ( NCYC, AC  12 ) )  CHA0066 

C  L AGEUL , A C 1 4 ) )  CHA0067 

C  NCYCPR, AC  17 ) )  CHA0068 

CSTAB,AC18)),CDTBA, A C 19 ) ) , C DTKOD, AC  20 ) ) , CKDT,AC21))  CHA006  9 

COMMON  /TOT/AMTO T . ETOT, EKTQT , FPTOT , TENT OT  CHA0070 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXX?*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA0071 
T=0 .  CHA00/2 


EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 


29 


-v. 


•  **- if- 


FILE i  CHARGEP  FORTRAN  AI 


NCYC«G 

JJJ*0 

CALL  NETUNM 
DELT*DT 
CALL  BEGIN 


1 

2 

X 

3 


CALL  SAFAE 


(L,X,U,P,R0,G,E,DU,DP,DR0,DG,DX5I,MIN, 

US,P5,UID0T,PID0T, 

FIMZ, ZMDOT  , 

TENA, FIRO, FIM, FIE, GIP, VOL, Z,DZ) 


1 

2 

x 

3 


CLpX,U,F,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

USp  PS, UI DOT , PI DOT , 

FIMZ, ZMDOT , 

TENA, FIRO, FIM, FIE, GIP, VOL, Z,DZ) 

1  NCYCaNCYC+l 
C  TIME  STEP  CONTROL. 

DT=DTBA 

IFCDT.GT . 1 . 1DOXDTKOD. AND. DTKOD.NE. 0 , )  DT*1 . 1DOXDTKOD 
IFINCYC . EQ . 2)  DT=DT/10.D0 
IF  (NCYC.EQ.l)  DT»0. 

IFCDT.EQ.Q.)  GO  TO  11 
NHAD=( (TMUD-T l/DT-l . D-10) 

IF(NHAD.GE.IO)  GO  TO  11 
DT*( TMUD-T )/DFLOAT  <  NHAD+1 ) 

CONTINUE 
T=T+DT 

I F( ( NCYC/ NCYCPR ) XNCYCPR . NE . NCYC . AND . NCYC . GT . NCYCPR)  GO  TO  33 
PRINT  10,  NCYC, T > DT , KDT 

FORMAT ( IX, *NCYC=*,I4,3X, ' T* • , Dll . 4, 3X, ' DT= » , Dll . 4, 3X, 'KDT"', 14) 
CONTINUE 
DTBA=DTMUD 
KDT  =  0 
NERI-1 

IF  ( DABS(T-TMUD) . LT . 1 . D-8 )  NERI=0 
CALL  CYCEUL 

L  ( L  ,  X,  U,  P ,  RO, G,  E,  DU,  DP,  DRO,  DG,  DXSI , MIN, 

!  US,PS,UIDQT, PI DOT , 

FIMZ, ZMDOT, 

TENA, FIRO , FIM, FIE, GIP, VOL , Z, DZ) 


11 


10 

33 


x 

3 


CALL  SAFAE 


(L,X,U,P,RO,G,E,DU,DP, DRO, DG, DXSI , MIN, 
US, PS, U I  DOT , PI DOT , 

FIMZ , ZMDOT , 

TENA, FIRO, FIM, FIE,GIP, VOL,Z, DZ) 

IF  CNERI.NE.C)  GO  TO  2 
CALL  PRINT 

1  (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

2  US, PS, UI DOT, PI DOT, 

x  FIMZ,ZMDOT, 

3  TENA, FIRO, FIM, FI E , GIP , VOL , Z , DZ) 

IF  ( DABS(T-TMUD) . LT . 1 . D-8 )  TMUD=TMUD+DTMUD 
CONTINUE 

DTKOD=DT 

IF  CT.LT.TMAX-1 .0-8)  GO  TO  1 
RETURN 
END 


Subroutine  NeTunM 

IMPLICIT  REAL*8(A-H,0-Z,$) 

COMMON  /AB/AC50) 

EQUIVALENCE  (L,A(1>) 

EQUIVALENCE  ( L L , A ( 2 ) ) , ( T , A ( 3 ) ) , ( DT , A ( 4 ) ) , ( TMAX , A 1 5 ) ) , 

1  (TMUD, A(6 ) J, (DTMUD, AC  7 ) ) , (JOB, A(8) ) , (NERI , A( 9) ) , 

2  (JJJ,A(10)), (KEYMON, A( 11 ) ) , (NCYC, A ( 12) ) 

EQUIVALENCE  ( COL ELA , A(13 ) ) 

EQUIVALENCE  ( L AGEUL , A( 1 4 ) ) 

EQUIVALENCE  ( KEYEK, A ( 16 ) ) 

EQUIVALENCE  ( NCYCPR , A (17 ) ) 

EQUIVALENCE  < STAB, A( 18 ) ) , C DTBA, A( 19 ) ) , ( DTKOD, A( 20 ) ) , ( KDT , AC 21 ) ) 
COMMON/DETO/QDET , PCJDET, RCJDET , UCJDET, DCJ DET , PODET , ROODET , 

1  RATE, TEMPC 

COMMON/DI FFUS/U2,PZ,R02,ARW 

COMMON  / DRAW/ GO DEL X, GO DEL Y, UMIN, UMAX, PMIN, PMAX, ROMIN, ROMAX 


CHA0073 

CHA0074 

CHA0075 

CHA0076 

CHA0077 

CHA0078 

CHA0079 

CHA0080 

CHA00S1 

CHA0082 

CHA0083 

CHA0084 

CHA0085 

CHA0086 

CHA0087 

CHA0088 

CHA0089 

CHA0090 

CHA0091 

CHA0092 

CHA0093 

CKA0094 

CHA0095 

CHA0096 

CHA0097 

CHA0098 

CHA0099 

CHA0100 

CHA0101 

CHA0102 

CHA0103 

CHA0104 

CHA0105 

CHA0106 

CHA0107 

CHA0108 

CHA0109 

CHA0110 

CHA0111 

CHA0112 

CHAC113 

CHA0114 

CHA0115 

CHA01 16 

CHA0117 

CHA0118 

CHA0I19 

CHA0120 

CHA0121 

CHA0122 

CHA0130 

CHAOiil 

CHA0132 

CHA0133 

CHA0134 


HSrvNtf 


CHA0136 
CHAO) 37 
CHAO] 38 
CHA0I39 
CHA0140 
CHAO  141 
CHA0142 
CHA0143 
CHA0144 
CHA  0145 
CHA0146 
CHA  0147 
CHA0148 
CHA0149 
CHA0150 
CHA0151 
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nnn  o  no 


CHAROEP  FORTRAN  A1 


1  »XMIN,XMAX, SMIN, SMAX, IVERSA 

COMMON  /GAM/GAMA, NO, MU2, 01, 02, G3,G4, 05, 06, 07, 08, 09, 010,011 

1  , 012,013, 014,015, 016,017, G1 8, G1 9, G20, 02 1,022,023 

2  , 024 , 025,026 , 027 , G28 ,029,030, 031 , 032,033,034 ,  G35 
REALMS  NO, MU2 

NAMEL 1ST  /IN/ LIN, GAMA , DT , TMUD , DTMU D , TMAX , 

1  OODELX, GODELY, UMIN , UMAX, PMIN, PMAX, ROMIN, ROMAX, 

2  SMIN, SMAX, IVERSA, KEYMON, COLEL A, STAB 

3  , LAGEUL,KEYEK 

4  , QDET 


CHA0152 
CHA0153 
CHA0154 
CHA0155 
CHA0156 
CHA0157 
CHAO  158 
CHA0159 
CHA0160 
CHA0161 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX3«XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA0162 


10 


LIN*L 

CHA0163 

LAGEUL-2 

CHA0164 

NCYCPR*1 

CHA0165 

KEYEK=1 

CHA0166 

TMUD=0. 

CHA0167 

DTMUD=1Q . DO 

CHA0168 

TMAX* 100, DO 

CHA0169 

STAB*0 . 5D0 

CHA0170 

DT-l.D-2 

CHA0171 

KEYM0N=1 

CHA0172 

GAMA=3 . DO+1 . D-6 

CHA0173 

QDET=C . 04D0 

CHA0174 

RATE*0 . 

CHA0175 

TEMPCSI . D50 

CHA0176 

GODELX-16DO 

CHA0177 

GODELY=20 . DO 

CHA0178 

IVERSAs100 

CHA0179 

UMIN*0 . 

CHA0189 

UMAX*  1 . DO 

CHA0181 

PMIN=0. 

CHA0182 

PKAX=0.5D0 

CHA0183 

ROMIN=0 . 

CHA0184 

ROMAX =3. DO 

CHA0185 

SMIN=0. 

CHA0186 

SMAX=0 . 03DO 

CHA0187 

COLELA=0. 

CHA0188 

READ  IN 

CHA0189 

CHA0190 

PRINT  IN 

CHAO  1 9 1 
CHAO 192 

GG=2 . DOXGAMA/ CGAMA-1 . DO ) 

CHA0193 

NG*GG 

CHA0194 

CONTINUE 

CHA0195 

MU2=(GAMA-1 . DO )/( GAMA+1 . DO ) 

CHA0196 

Gl=GAMA~i . DO 

CHA0197 

G2  =  l . D0-MU2 

CHA0198 

G3=2.D0/(3. DOXGAMA-1 . DO ) 

CHA0199 

G4*( GAMA+1 . DO )/2 . DO 

CHA0200 

G5  =  0 .5D0*( 3. DOXGAMA-1 .DO)/ (GAMA+1 .DO) 

CHA0201 

G6  =  ( GAMA+1 . DO )/(2 . DOXGAMA) 

CHA0202 

G7=2 . DO/ (GAMA- 1 .  DO ) 

CHA0203 

G8  =  (GAMA-.l  .D0)/(2. DOXGAMa) 

CHA0204 

G9=(GAMA+1 .D0)/(4. DOXGAMA) 

CHA0205 

G10=l . DO/GAMA 

CHAQZ06 

Gll=( GAMA+1. DQ)/4. DO 

CHA0207 

G12=GAMA/( GAMA-1 . DO) 

CHA0208 

G13  =  0.5D0X( GAMA- 3  -  DO )/( GAMA+1 . DO ) 

CHA0209 

G14=0.5D0X(3. DOXGAMA-5. DO)/ (GAMA+1 .DO) 

CHA0210 

G15=GAMAX( 3. DOXGAMA-1 .DO) 

CHA0211 

G16=( GAMA+1 . D0)/<2 . DO X( GAMA-1 . DO )  ) 

CHA0212 

G17  =GAMA+1 . DO 

CHA0213 

G18=GAMA*( GAMA+1 .DO)/ (3. D0*GAMA-1 -DO) 

CHA0214 

G19  =  ( 3. DOXGAMA-1 . DO)/ (GAMA+1  .  DO  ) 

CHA0215 

G2Q=2 . DOX( GAMA-i . DO )/( 3 . DOXGAMA-1 . D0)**2 

CHA02I6 

G?1=GAMA*( 3. DOXGAMA-5. DO )/( 3 . DOXGAMA-1 . D0)X*2  CHA0217 

GODELX=GODELX/2.54DO 

CHA0213 

GODELY* GODELY/2. 54 DO 

CHA0219 

CALL  NAMPL T ( IVERSA ) 

CHA0220 

CALL  LIMIT(IOOO.DO) 

CHA0221 

CALL  PLOT( 0 . , 0 . 5D0 , -3 ) 

CHA0222 

P0DET  =  0  . 

CHA0223 

FILE.  CHARGEP  FORTRAN  A1 


I.' 


ROODET*0. 

PCJDET-0. 

UCJDET»0. 

DCJDET-O. 

RCJDET*0. 

IFCQDET . LE. 0 . )  GO  TO  100 
C  DETONATION  DATA 
0LET*0 . Q4D0 
P0DET«0. 

R00DETS1 .8D0 

PCJDET=PODET-(GAMA-1.DO)*(-QD£T)*ROOD£T+ 

1  DSQRTC<(GAMA-1.D0)XQDETXRQQDET)XX2-2.D0*MU2XGAMAX 

2  (-QDET)XPODETXROODET) 

RCJDET‘RC0DETX((GAMA+1.DG)XPCJDET-P0DET)/(GAMAXPCJDET) 

CCJ  »DSQRTC  GAMAXPCJ  DET/RCJ DET ) 

DCJDET=CCJXRCJDET/ROODET 
>JCJDCT*DCJDET~CCJ 
PRINT  101 

101  F0RMATC1H1,/, IX, 'DETONATION  DATA'/) 

PRINT  102,  QDET, GAMA, TEMPO, RATE 

102  FORMAT (/IX, ’ QDET,GAMA, TEMP, RATE= • , 4D18 .8) 

PRINT  103,  ROODET , PODET 

103  F0RMATC/1X, 'UNBURNED  STATE  ROODET, PODET** ,2D18 .8) 

PRINT  104,  PCJDET, PCJDET ,RCJDET ,UCJDET 

104  F0RMATC/1X,  'CJ  POINT  DCJDET, PCJDET, RCJDET, UCJDET= • , 4D18 . 8) 

100  CONTINUE 

RETURN 


CHA0224 

CHA0225 

CHA0226 

CHA0227 

CHA0228 

CHA0229 

CHA0230 

CHA0231 

CHA0232 

CHA0233 

CHA0234 

CHA0235 

CHA0236 

CHA0237 

CHA0238 

CHA0239 

CHA0240 

CHA0241 

CHA0242 

CHA0243 

CHA0244 

CHA0Z45 

CHA0246 

CHA0247 

CHA0248 

CHA0249 

CHA0250 


1  ( I. ,  X,  U,  P ,  RO , G,  E,  DU,  DP,  DRO,  DG,  DXSI  ,MIN,  &E6M 

2  US,  PS, 'J I  DO  V,  PI  DOT, 

X  FIMZ,ZMDOT , 

3  TENA, FIRO,FIM,FIE,GIP,VOL,Z,DZ) 

IMPLICIT  REALX8(A-H,0-Z,$) 

DIMENSION  X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L), DRO( L ) , 


CHA0252 
CHA0253 
CHA0254 
CHA0255 
CHA0256 
CHA0257 
CHA0258 

1  DGC  L ) ,  DXSK L) , MINC L ) ,  CHA0259 

2  US( L ) , PS( L ) , UIDOT ( L ) , PIDOT ( L )  CHA0260 

3  , TENAC L ) , FIROC L ) , FIMC t ),FIE(L)  CHA026I 

4  ,GIPCL),VOL(L),Z(L),DZ(L)  CHA0262 

5  , FIMZ(L) ,ZMDOT(L)  CHA0263 

COMMON  /AB/AC50)  CHA0264 

COMMON  /GAM/GAMA , NG, MU2, G1 , G2,G3 , G4 , G5, G6 , G7 , G8 ,G9, G10 ,GI1  CHA026 5 

1  ,GI2,G13,GI4,G15,G16,G17,G18,GI9,G20,G21, G22 . G23  CHA0266 

2  , G24 , G25 , G26 . G27 , G28 , G29 , G30 , G31 , G32 , G33 , G34 , G35  CHA0267 

REALX8  NG , MU2  CHA0268 

EQUIVALENCE  (LL,A(2))  CHA0269 

EQUIVALENCE  < LAGEUL , AC  14) )  CHA0270 

EQUIVALENCE  (UGAL,A(15))  CHA0271 

EQUIVALENCE  (STAB,  AC  18) ) ,  CDTBA,AU9)) ,  ( DTKOD,  AC 20 ) ) ,  C KDT ,  AC21 ) )  CHA0272 

COMMON/DETO/QDET, PCJDET, RCJDET, UCJ DET, DCJDET, PODET, ROODET,  CHA0273 

i  RATE, TEMPO  CHA0274 

COMMON  / DRAW/GO DEI X, GO. )ELY, UMIN, UMAX, PMIN, PMAX, ROMIN, ROMAX  CHA027  5 

1  ,XMIN,XMAX,SMIN,SMAX, IVERSA  CHA0276 

COMMOM/GIT/ROLIM, 3L IM, XGITC 200 ) , ROGIT ( 200 ) , ROUGITC 200 ) ,EGITC200)  CHA0277 
COMMON  /GITN/NPO  CHA0278 

LOGICAL  CSOF  CHA0279 

CXXXi«XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX*XXXXXXCHA0280 
DT3A=0 .  CHA0281 

DTKOD=0,  CHA0282 

KDT  =  0  CHA0283 

P0=l.D-9  CHA0284 

RHOO=1 . D-7  CHA0285 


UO  =  UCJ  DET 
UGAL  =0 . 

X0  =  0  . 

XI =50 . DO 
XCHARG= 1 0 , 
XMIN=XO 
XMAX=X1 


CHA0286 

CHA0287 

CHA0288 

CHA0289 

CHA0290 

CHA0291 

CHA0292 


CHA0293 

CHA0294 

CHA0295 


v  v  s."  •/  •  /  ■  *  •  •  s*  -  *  •  ■ 


ooo 


CHARGEP  FORTRAN  A1 


31 


32 

3Q 

2 


CONTINUE 

X(L)»X1 

U0«U0*CR0SS(XCHARG) 

DO  2  1*2, LL 

IFCI.GT.2)  U( I ) *U0/CR0SS(X( I ) ) 

P(I)*PQ 

RO(I)*RHOO 

Z(I>*0. 

00  TO  (31,32),  LAGEUL 
CONTINUE 

E(I>*PCI)/<(GAMA-1.DO>*RO<I))+0.5DO«U(I)*X2+Z(I)*QDET 

GO  TO  30 

CONTINUE 

E(I>«F(I)/(GAMA-1.DO)+0.5DO*RO<I)XU<I)XX2+ZCI)*RO(I)*QDET 

CONTINUE 

G( I ) *  DSQRT  < GAMAXp ( I ) XRO ( I ) ) 

CONTINUE 

DO  3  1*2, LL 

TENA(I)=RO(I)XU(I) 

VOL(I)*(X(I+1)-XCI))X(X(I+1)XX2+X(I+1)XXCI)+X(I)X*2)/3.DO 

CONTINUE 


INSERT  DETONATED  CHARGE  FLOW  FIELD  FROM  TAYLOR'S  SOLUTION. 


CALL  TAYLOR (GAMA) 

RONORM*RCJDET 
RUNORM*RCJDETXUCJDET 
ENORM=RCJ DETXDCJ  DETXX2 
XLIM=XGIT(NPO) 

NGIT=NP0-1 

XG1*XLIM 

XG2=XGIT(NGIT) 

AROIP  =ROGIT  (NPO)+ROLIMXXLIMXX3/3.DO 
AROUIP=ROUGIT(NPQ) 

AEIP  =EGIT  (NPO)+  ELIMXXLIMXX3/3.D0 

XP=X(2)/XCHARG 

DO  100  1=2, LL 

IP=I+1 

XI=XP 

AROI  =AROIP 
AROUI=AROUIP 
AEI  =AEIP 
XP=X( IP  J/XCHARG 

IF( DABS( XP-1 . D0).LT.1.D-10)  XP=1.D0 
CS0F=(XP.GE.1.D0) 

IF(XP.GE.XLIM)  GO  TO  101 
:  UNIFORM  FLOW  REGION 

DELV0L=(XLIM-XP)*(XLIMxx2+XLIMXXP+XPXX2)/3. DO 
AROIP  =ROGIT  (NPO)+ROLIM*DELVOL 
AROUIP=ROUGIT ( NPO ) 

AEIP  =EGIT  ( NPO )+  ELIMXDELVOL 
GO  TO  102 
101  CONTINUE 


104 


103 


NON  UNIFORM  FLOW  REGION. 

I F( . NOT . CSOF)  GO  TO  104 

LAST  POINT.  (THIS  IS  THE  DETONATION  FRONT  POINT 
AROIP=  0. 

AR0UIP  =  0  . 

AEIP=  0. 

GO  TO  102 
CONTINUE 

I F(XP . LE . XG2)  GO  TO  103 
NGIT  =  NGI T-l 

IF(NGIT.LE.O)  CALL  SOF( ' BEGIN  104.  NGIT.LE.O.') 
XG1 =XG2 

XG2=XGIT ( NGIT ) 

GO  TO  104 
CONTINUE 

FRAC=(XP-XG1)/(XG2-XG1) 

IF(FRAC. LT . 0 .  )  CALL  SOF( 'BEGIN  103. 

IF(FRAC.GT.l.DO)  CALL  SOF( 'BEGIN  103. 


X  =  l)  . 


FRAC.LT.O.  '  ) 
FRAC.GT.l.  ') 


AROIP  =( 1 . DO-FRACJXROGIT  ( NGIT+1 )+FRAC*ROGIT  (NGIT) 


33 


CHAQ296 

CHA0297 

CHA0298 

CHA0299 

CHA0300 

CHA0301 

CHA0302 

CHA0303 

CHA0304 

CHA0305 

CHAC306 

CHA0307 

CHA0308 

CHA0309 

CHA0310 

CHA0311 

CHA0312 

CHA0313 

CHA0314 

CHA03J5 

CHA0316 

CHA0317 

CHA0318 

CHA0319 

CHA0320 

CHA0321 

CHA0322 

CHAO 32 3 

CHA0324 

CHA0325 

CHA0328 

CHA0327 

CHA0328 

CHA0329 

CHA0330 

CHA0331 

CHA0332 

CHA0333 

CHA0334 

CHA0335 

CHA0336 

CHA0337 

CHA0338 

CHA0339 

CHAQ340 

CHA0341 

CHA0342 

CHA0343 

CHA0344 

CHA0345 

CHA0346 

CHA0347 

CHA0348 

CHA0349 

CHA0350 

CHA0351 

CHA0352 

CHA0353 

CHA0354 

CHA0355 

CHA0356 

CHA0357 

CHA0358 

CHA0359 

CHA0360 

CHA0361 

CHA0362 

CHA0363 

CHA0364 

CHA0365 

CHA0366 

CHA0367 


ooo 


FILE)  CHAROEP 


FORTRAN  A1 


D0-FRAC)XROUO2TCNGIT+l)<tFRACHR0UGITCNOlT) 
DO-FRAOXEOIT  !NQIT+1)+FRACXEGIT  (NOIT) 


MOMENTUM 

DO 


AR0UXP*!1 

aeip  «u 

102  CONTINUE 

C  COMPUTE  MASS*  MOMENTUM  AND  ENERGY  DENSITIES. 

IFIXP.LE.XLIM)  GO  TO  105 
C  CONSERVATION-FORM  DEFINITION  OF  MASS* 
DV0L><XP-XI)x(XP*X2+XPxXI+XImk2)/S 
RO  (I)"RONORMX(AROI  -  AROIPVDVOL 
TENA(I)-RUNORMX(AROUI-AROUIP)/DVOL 
E  (D-ENORM  *CAEI  -  AEIP)/DV0L 
GO  TO  106 
105  CONTINUE 
C  UNIFORM  FLOH  REGION 

RO  ( I ) *R0NQRM*R0L IM 
TENA! I)*0 . 

E  (D-ENORM  X  ELIM 
CONTINUE 

U!I)-TENA!I)/RO(I) 

P(I)-(GAMA-1.DO)X(E(I)-0.5DOXRO(I)XU(I)XX2) 
PRINT  m,I,CSOF,U(I),P(I),RO(I),E<I) 
FORMAT!/ IX* * I ,CSOF, U, P,  RO, E* ' ,  14, L3, 4D14 .4) 
IF(CSOF)  GO  TO  109 
CONTINUE 
CONTINUE 
DO  A  1*2, LL 

PXSI<I)*<X(I+l)-X(I))XRO<I) 

CONTINUE 
RETURN 


AND  ENERGY  DENSITY. 


106 


111 


100 

109 


"SffcniDTnjr 


TAJ  L  Off 


CHA0368 
CHA0369 
CHA0370 
CHA0371 
CHA0372 
CHAC373 
CHA0374 
CHA0375 
CHA0376 
CHA0377 
CHA0378 
CHA0379 
CHAO 380 
CHA0381 
CHA0382 
CHA0383 
CHA0384 
CHA0385 
CHA0386 
CHA0387 
CHA0385 
CHA0389 
CHA0390 
CHA0391 
CHA0392 
CHA0393 
CHA0394 
CHA0395 
CHA0396 


TAYLflftfSSWAl - 

INPLICIT  REALMS! A-H, O-Z, ♦ ) 


TAYLOR  —  SELF  SIMILAR  SPHERICAL  DETONATION  !CJ)  FLOH  FIELD 


CHA0397 
CHA0398 
CHA0399 
CHA0400 
CHA0401 

COMMON  /GGGG/G,G1,G2,G3,G4,G5,G6,G7,G8, G9,G10  CHA0402 

COMMON  /PAR/RHOO , QO , ROC J , DCJ , UCJ , PC J , DPSI , PSIMAX, CO , UO  CHA0403 

COMMON  /GITN/NPO  CHA0404 

COMMON/GIT/ROLIM,  ELIM,  XGITC200  ) ,  ROGIT(200 ) ,  ROl'GIT!  200  ),EGIT!  200)  CHA0405 

CXXXXXKKXXXKXXKXXXXXXXXX*XSXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXKXXXXXXXXX*MXXCHA0406 


G*GAMA 

CHA0407 

PRINT  101 

CHA0408 

101 

FORMAT! *1') 

CHA0409 

PRINT  110 

CHA0410 

110 

FORMAT ( IX, 'G .  I.  TAYLOR  SOLUTION.  N,PSI,U,C,X/AM,AT,AE3'//) 

CHA0411 

CALL  INIDAT 

CHA0412 

X-l . DO 

CHA0415 

Y-0 . 

CHA0414 

U*UO 

CHA0415 

C3CO 

CHA0416 

AM=0. 

CHA0417 

AT-O. 

CHA0418 

AE  =  Q. 

CHA0419 

PSI*-DLOG(U) 

CHA0420 

DO  1  N3 1 , NPO 

CHA0421 

XGIT  ! N  )  3X 

CHA0422 

R0G1T  (N)=AM 

CHA0423 

ROUGH!  N )  =AT 

CHA0424 

EGIT  (N)=AE 

CHA0425 

PRINT  11,  N,PSI,U,C,X,AM,AT,AE 

CHA0426 

il 

FORMAT! 1X,I4,4D14.5/5X,3D14.5) 

CHA0427 

CALL  RUNGE!N,PSI,X,C,AM, AT , AE, PSIN , XN, CN, AMN, ATN, AEN) 

CHA0428 

PSI =PSIN 

CHA0429 

U=DEXP< -PSI ) 

CHA0430 

X  =  XN 

CHA0431 

C=CN 

CHA0432 

AM  =  AMN 

CHA0433 

AT “ATN 

CHA0434 

AE=AEN 

CHA0435 

1 

CONTINUE 

CHA0436 

RQLIM=!C/C0HXG3 

CHA0437 

ELIM=G5X!C/C0)X*G4 

CHA0438 

AM03AM+!C/C0)xxG3xX»-<3/3.D0 

CHA0439 

34 
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AM0»AMO*5 . DO*(G*l . DO )/G  CHA0440 

AE0“AE+(GIixCC/C0)XXG4+0. 5D0X<C/CO)XXCS3XUXX2)XXXX3/5.DO  CHA0441 

AEO»AEOX6.DOX<Gn,DO)X(G»*2-\.DO)/G  CHAQ442 

PRINT  22,AM0,A£O  CHA0443 

22  F0RMAT(///1X, 'MASS  AND  ENERGY  CHECK  (SHOULD  BE  I.)*//  CHA0444 

1  IX, 'MO3 ' , D17 . 8 , SX, ' E0“ ' , D17 . 8// )  CHA0445 

RETURN  CHAQ446 

- iuSaCUTINETNIBST -  - m*r - MM 

IMPLICIT  REAL*3(A-H,0~Z,*>  CHA0449 

COMMON  /GGGG/G, G1,G2, 83 ,G4,G5,G6, G7,G8,G9,G10  CHAC450 

COMMON  /PAR/RHQO , QG , ROCJ , DC J , UCJ , PCJ , DPSI , PSIMAX, CO , UO  CHA0451 

COMMON  /GITN/NPO  CHA0452 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA0453 
NP0*20Q  CHA0454 

PSIMAX*IO . DO  CHA0455 

UO«1.DO/(G+1.DO)  CHA0456 

C0»1.D0-U0  CHA0457 

DP5I*PSIMAX/DFL0AT(NP0)  CHA0458 

G1*G~1 . DO  CHA0459 

G2*G1/2.DP  CHA0460 

G5*2. DO/CG-l . DO)  CHA0461 

G4«2.D0*G/CG-1.DQ)  CHA0462 

G5»G/<(G+1.DO)XX2X(G-1.DO))  CHA0463 

RETURN  r  CHA0464 


^fMOmrTnjfiSETTrnSSITjC;  C ,  A  M ,  A  T ,  A  E ,  P  3 1 N ,  X  N ,  C  N ,  A ,  1 N ,  A  T  N ,  A  E  N  ) 
IMPLICIT  REALX8(A-H,0-Z,*) 

COMMON  /GGGG/G , G1 , G2 , G3 , G4 , 05 , G6 , 07 , G8 , G9 , G1 0 
COMMON  /PAR^RHOO , QO , ROCJ , DCJ , UCJ , PCJ , DPSI , PSIMAX, CO , UO 
COMMON  /GITN/NPO 


CHA0466 

CHA0467 

CKA0468 

CHA0469 

CHA0470 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA0471 


H=DPSI 
H2=H/2 . DO 
H6=H/6 . DO 

CALL  DERIV(PSI,X,C,AM,AT,AE, 

1  DXDP1 , DCDP1 , DMDP1 , DTDP1 , DEDP1 3 

CALL  DERIV(P3I+H2,X+H2KDXDPl,Ci-H2*DCDPI,AM,AT,AE, 

I  DXDP2, DCDP2, DHPP2, DTDP2, DEDP2) 

CALL  DERIV(PSI+H2,X+H2XDXDP2,C+-H2*DCDP2,AM,  AT,AE, 

1  DXDP3, DCDP3, DMDP3, DTDP3,DEDP33 

CALL  DERIVC  PSI+H, X+H*DXDP3, C+HXDCDP3 , AM , AT, AE, 

1  DXDP4, DC DP 4, DMDP4, DTDP4, DEDP4 ) 

PSIN=PSI+H 

XN  =  X+H6*( DXDP1+2 . DQX( DXDP2+DXDP3  3  +  DXDP4  3 

CN=C+H6X( DCDP1  +  2 . DO*U  DCDP2+DCDP3  3  +  DCDP4  3 

AMN=AM+H6*(DMDP1+2.D0*(DMDP2+DMDP33+DMDP43 

ATN=AT+H6*(DTDP1+2.D0X(DTDP2+DTDP33+DTDP43 

AEN=AE+H6X( DEDP1+2 .D0X(DEDP2+DEDP33+DEDP43 

RETURN  ~ 

end  <y 

SUBROUTINE  Dt'RTVTPSI ,  X ,  C ,'ffl7ST7AE7DXD P 7DCDP ,  DMDP ,  DTDP ,  DEDP 3 
IMPLICIT  REAL#8(A-H,0-Z,$) 

COMMON  /GGGG/G, G1,G2.G3,G4,G5,G6,G7,G8,G9,G10 
COMMON  /PAR/RHOO,QO, ROCJ, DCJ, UCJ, PCJ, DPSI , PSIMAX, CO , UO 
COMMON  /GITN/NPO 


CHA0472 
CHA0473 
CHA0474 
CHA0475 
CHA0476 
CHA0477 
CHA0478 
CHA0479 
CHAO 480 
CHAO 481 
CHA0482 
CHA0483 
CHA0484 
CHA0485 
CHA0486 
CHA0487 
CHA0488 
CHA0489 
_CHA.QA2.fl — 
CHA0491 
CHA0492 
CHA0493 
CHA0494 
CHA0495 


CX>:>XXXXXXXXXXX?<XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX**XX*XXXXXXXXXXXXXXXX*X*XCHA0496 


U=DEXP( -PSI 3 

DXDP=0.5DQ*Xx(C-U+X3*(C+U-X3/CXX2 

DCDP=-G2XUX(X-U3/C 

DMDP  =  -(C/CO  3XXG3XXXX2XDXDP 

DTDP=DMDP*U 

DEi3P  =  -rG5*(C/C0)**G4+0 . 5D0XCC/C0 3X*G3XUX*23*XXX2*DXDP 
RETURN 


DOUBLE  PRECISION  FUNCTION  RATI0CX3  RATIO  CHA0505 

IMPLICIT  REAL*8(A-H,0-Z,$)  CHA0506 

CXXXXX#X**XXXXXXXXXXXXXX;*##X*XXXXXXX*XXXXXXXXXX#XXXXXX#XXXXXX#XXXXXXXXXXCHA0507 
RATI0=0.  CHA0508 

IF(X.LE.l .D-83RETURN  CHA0509 

RATIQ=2 . DO/X  CHA0510 

RFTURN  CHA0511 


CHA0497 

CHA0498 

CHA0499 

CHAO5O0 

CHAQ501 

CHAC502 

CHA0503 

CHA0SQ4- 

CHA0505 

CHA0506 


r 
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I 
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■£tUL 


DOUBLE  PRECISION  FUNCTION  CROSS(X)  CROSS  CHA05I3 

IMPLICIT  REALX8(A-H,0-Z,$>  CHA0514 

CXKXXXXXXXXXXXXXXXKKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXCHA0515 


CROSS'l . DO 
CR05S=XXX2 
RETURN 


END _ 

SUfWOTTWF 


cjCrOiS 


CWfeUL 

1  (  L ,  X,  U,  P, RO, 0, E, DU, DP, DRO, DO/ DXSI ,  MIN, 

2  US, PS, UIDOT , PI DOT , 

X  FIMZ, ZMDOT , 

3  TENA, FIRO, FIM, FIE, GIP, VOL ,Z, DZ) 

IMPLICIT  REALX8L A-H, 0-Z, $ ) 

DIMENSION  X(L),U(L),P(L),ROCL),G(L),E(L),DU(L),DP(L), DROCL ) , 

1  DGC  L ) , DXSI ( L ) , MINC  L ) , 

2  USCL),PS(L) , UIDOT  <  L ) , PIDQT( L ) 

3  ,TENA(L),FIRO(L),FIM<L),FIECL) 

4  , GIP( L ) , VOL (L),Z(L), DZC L ) 

5  ,fimz(d,zmdot:l) 

COMMON  /AB/A(50) 

EQUIVALENCE  ( LL , A(2> ) , (T, A( 3 ) ) , ( DT, AC  4) ) , CCOLELA, AC  1 3) ) 
EQUIVALENCE  ( KEYEK, AC  16 ) ) 

EQUIVALENCE  (STAB, AC18 ) ) , ( DTBA, A( 19 ) ) , ( DTKOD, A(20 ) ) , ( KDT, A(21 ) ) 
COMMON  /GAM/GAMA, NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10, Gil 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2  ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35 
REAL *8  NG,MU2 

COMMON  /TOT/AMTOT, cTOT, EKTOT, EPTOT, TENTOT 
COMMON  /AZOV/ I SAFA, NORIMN, USAF , PSAF, ROSAF, GSAF, ESAF, DPSAF 
1  , DXSIL , DXSIR 

LOGICAL  NORIMN 

COMMON  /STEPO/UL,PL»ROL,GL,UR,PR, ROR, GR, USTAR, PSTAR, 

1  RSTARL , RSTARR, GSTARL , GSTARR, 

2  CL , CR, CSTARL , CSTARR , SL , SR , WL , WR, UW( 6 ) 

3  , LAMDAL , LAMDAR, RATEL , RATER, TEMPL , TEMPR, TEMPSL , TEMPSR 

4  , ZL , ZR, ZSTARL , ZSTARR, NFLUX, HELEML , HELEMR 
REALX8  LAMDAL, LAMDAR 

LOGICAL  HELEML, HELEMR 

COMMON  /STEP1/DUIDT, DPIDT, DGIDTL , DGIDTR, DRIDTL  ,  DRIDTR 

2  , ASTARL , ASTARR, LAMDSL , LAMDSR, DSDAL , DSDAR, DZDAL , DZDAR 

3  , RAT , SH 

4  , BETACL , BETACR , DSDASL , DSDASR , DZDASL , DZDASR 
REALX8  LAMDSL , LAMDSR , DSDAL , DSDAR , DZDAL , DZDAR 

COMMON  /GRADS/ DUDXIL, DPDXI L, DGDXI L, DRDXI L , DZDXI L , DSDXI L , 

1  DUDXIR, DPDXI R, DGDXIR, DRDXIR , DZDXIR, DSDXIR 

COMMON  /FI/FIHI,FIH2,FIH3,UXN,PXN,GXN,R0XN,ZXN 

1  ,  GIH 

2  , FIH4 , ZMDOTL , ZMDOTR 

COMMON/ DETO/QDET , PCJDET , RCJDET , UCJDET , DCJDET, PODET, ROODET, 

1  RATE, TEMPO 

COMMON/ PULS/ PR ESS ( 10),PULSEI(10),PULSE2(10),PULSE3(10),PULSE4(10) 
DATA  ERRP/O.DO/ 

DATA  NERRP/O/ 

DATA  K0TZ/7777777777B/ 


CHA0516 
CHA0517 
CHAO  5  18 
GHAtlSiq 


DT2=DT/2 . DO 

CHA0568 

UXN  =  0  . 

CHA0573 

PXN  =  0 . 

CHA0574 

1 

ROXN=0. 

CHA0575 

j 

ZXM  =  0  . 

CHA0576 

DO  1  1=2, L 

CHA0577 

: 

H 

IM=I-1 

CHAO 57 8 

A 

UXNM=UXN 

CHA0579 

* 

PXNM=PXN 

CHA0580 

X 

ROXNM=ROXN 

CHA0581 

ZXNM=ZXN 

CHA0582 

i 

Ul=U(IMM0.5D0*DU(IM) 

CHA0583 

PL=P(IM)+0.5D0*DP(IM) 

CHA0584 

ROL=RO(IM)+0.5DO*DRO(IM) 

CHAU585 

GL  =DSQRT (GAMA*PL*RQL ) 

CHA0586 

* 

CL=GL/ROL 

CHA0587 

mmamm 


CHARGEP  FORTRAN  A1 


PAGfJ 


ZL*Z11M)+0.,?DO*D2<IM) 

IFCZL.GT.l.DG)  £1*1 . DU 
IFC2L.LT. 0.  )  ZL»0 

SL-PL/CQIXROLXXOAKA) 

TEMPI «PL/R0L 
UR*U(I)-0.300*DU(I) 

PR*P(I)-0.5D0XDP(I) 

ROR»ROa)~0.5DOKORO(I) 

GR- OSQRT ( GAMAXPRXROR ) 

CR-GR/ROR 

ZR*Za)“0.5D0XDZJI) 

IFCZR.GT.1.D0)  ZR*1 .DO 
IFCZR.LT. 0.  )  ZR*0. 

SR“PR/ C  G1 XRORXNGAMA ) 

TEMPR*PR/ROR 

C 

CALL  KXEMAN(L,I,MIN) 

C 

DUDXIL*DUCIM)/DXS1(IM) 

DPDXIL*DP(IM)/DXSICIM) 

DRDXIL»DRQ(1M)/0XSI(1M) 
DGDXIL«0.5!)0XGLX(DPDXIL/PL+0RDXIL/R0L) 
DZDXIL“DZ(IM)/DXSIC IM) 
DSDXIL*SLX(DPDX1L/PL-GAMAXDRDXIL/R0L) 
DUDXIR*DJ( I )/DXSI ( x ) 

DPDXIR*DPCI)/DXSICI) 

DRDXIR»DRO(I J/DXSI Cl ) 
DGDXIR*0.5D0XGR»iCDPDXIR/PR+DRDXIR/ROR) 
DZDXIR*DZ(I  VDXSICI) 

DS0X1R*SRX( DPDX1R/PR-3AMAXDRDX1R/ROR) 
SHKCROSS(XC I) ) 

RAT*RATIO(X(I)J 

C 

CALL  MAGA(L,I,MIN) 

C 

USC  2 )*USTAR 
PS( I )*PSTAR 
UIDOT(I)=DUIDT 
PID0T(I)*DPID7 
C 

CALL  FLUXEC L , I , MIN) 

C 

FXROC I )=FIH1 
FIM  (I)=FIH2 
FIE  ( I ) SFIH3 
FIMZ( I ) =FIH4 
GIPC I ) =GIH 
DUCIM)=UXN-UXNM 
DPC IM) =PXN-PXNM 
DRO(IM)=ROXN-ROXNM 
DZ(IM)=ZXN-ZXNM 
C  STATIONS  OUTPUT 

IF((I-42)*(I-62)*(I-82mi-102)  .NE.Cl)  GO  TO  1 
NPU  =  0 

IF(I.EQ.42>  N?U  =  1 
IF(I . EQ.62)  NPU=2 
IFCI.EQ.82)  NPU=3 
IF(I.EQ.102)NPU=4 

IFCNPU.EQ.O)  CALL  SCFC'FLUXE  90.  NPU.EQ.O*) 
PRESSC  NPU )=GIH+FIH2 
PULSE1(NPU)=PULSE1(NPU)+DT*GIH 
PULSE2(NPU)=PULSE2(NPU)+DT*<GIH+FIH2) 
PULSE3(NPU)=PULSE3(NPU)  +  DT*FIH1*CR0SS(XCI) ) 
PULSE4(NPU)=PULSE4(NPU)+DTXFIH2*CRQS5(X( I) ) 

.  1  CONTINUE 
C 

AMTOT=Q. 

ETOT=0 . 

EKTQT=0 . 

EPTOT=0. 

TENTOT  =  0 . 

FI1=FIR0C2) 


CHA0588 

CHA0389 

CHA0590 

CHA9591 

CHA0592 

CHA0593 

CHA0594 

CHA0595 

CHA059$ 

CHA0597 

CHA0598 

CHA0599 

CHAOGOQ 

CHA0601 

CHA0602 

CHA0603 

CHA0604 

CHA0605 

CHA0606 

CHA06G7 

CHA0608 

CHA0609 

CHA0610 

CHA0611 

CNA0612 

CHA0613 

CMA0614 

CHA061S 

CHA0616 

CMA0617 

CHAC618 

CHA0619 

CHA0620 

CHA0621 

CHA0622 

CHA0623 

CHA0624 

CHA0625 

CHA0626 

CHA0627 

CHA0628 

CHA0629 

CHA0630 

CHA06  31 

CHA0632 

CHA0633 

CHA0634 

CHAQ635 

CHA0636 

CHA0637 

CHA0638 

CHA0639 

CHA0640 

CHA0641 

CHA0642 

CHA0643 

CHA0644 

CHA0645 

CHA0646 

CKA0647 

CHA0648 

CHA0649 

CKA0650 

CHA0651 

CHA0652 

CHA0653 

CHA0654 

CHA0655 

CHAC656 

CHA0657 

CHA0658 

CHA0659 
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C 

C 

291 

C 

29 

2 

C 

200 


FI2-FIN  (2) 

FX3*FIE  t 2) 

FX4*FIMZ<2) 

OI2«OIP(2) 

SH*CR0S5( X(2) ) 

DO  2  1*2*  f.L 

IP*I+1 

FIM1BFI1 

FIM2*F12 

FIM3»FI3 

FIM4»FI4 

GIM2»GI2 

SHM-SH 

FI1*FIR0( IP) 

FI2*FIM  (IP) 

FIS*FIE  (IP) 

FI4*FIMZ(IP) 

GI2-GIP  (IP) 

5H*CRQSSCX(IP)) 

DVOL*VOL( I ) 

ROOLD-RO(Z) 

POLD*P(I) 

EOLD*E(I) 

UOLD*U(I) 

Z0LD*Z(I) 

ZKODM*ZOLD*ROOLD 

T01D»P0LD/RQ0LD 

DX*X(IP)“X(I) 

DTVOL-DT/DVOI. 

RO(I)»RO(I)-DTVOL*(SH*FI1-SHM*FIM1) 

TENA(I)BTENA(I)-DTV0L*(SH*FI2-SHM*FIM2)-(DT/DX)X(GI2-GIM2) 

E(I)BE(I)-DTV0L*(SH*FI3-SHM*FIM3) 

U(I)BTENA(IVRO(I) 

Z( I ) B(ZKODM-DTVOL*( SH*FI4-SHM*FIM4) )/RO( I ) 

IF(Z( I ) .GT , 1 . DO )  Z(I)*1.D0 
IF(Z(i).LT.O.)  Z(I)*0. 

UAV=U( I ) 

ROAV=RO( I ) 

EP=E(I)-Q . 5D0#R0AV#UAV*#2 

IF(EP . GT . 0 . )  GO  TO  291 

NERRP=NERRP+1 

ERRP=ERRP+(l.D-3~EP)*DV0l 

IFC  ERRP . GT . 0 . 24U0)  GO  TO  291 

EPS1 . D-8 

CONTINUE 

IF(EP.LE.O.)  GO  TO  7001 
P(I>*G1#IE.P 

G( I ) =DSQRT(GAMA*P< I ) *R0( I )  ) 

UPC=DABS(U(I))+G(I)/RO(I) 

DTI--<>7AB*DX/UPC 

IF( DTI . GT . DTBA)  GO  TO  29 

DTBA=DTI 

KDT  =  1 

CONTINUE 

DXSI(I)*RO(I)*DX 
ETOT -ETOT+E( I )*DVOL 
EPTOT=EPTOT+EP*DVGL 
AMTOT =AMTOT+RO( I )*DVOL 
TENTCT=TEN10T+TENA( I )*DVOL 
CONTINUE 

EKTOT  =  ETOT-EPTO  T 

IFCCOLELA . cQ . 0 . )  GO  TO  200 
CALL  DCUL  E( L  *  X, U  , DU  ,MIN,1) 

CALL  DC0LE(L,X,P,DP,MIN,2) 

CALL  DCULcLL,X.R0,DR0,ttIN,3) 

CALL  DCCLE(L,X,Z,DZ,MIN,4) 

CONTINUE 

CALL  BD0X1 (L* X- U  , DU  ,MIN*1) 


CHA0660 
CHA0661 
CHA0662 
CHA0665 
CHA0664 
CHA0665 
CHA0666 
CHA0667 
CHA0S68 
CHA0669 
CHA0670 
CHA3671 
CHA0672 
CHA0673 
CHA0674 
CHA0675 
CHA0676 
CHA0677 
CHA0678 
CHA0679 
CHA0680 
CHA0681 
CHA0682 
CHA0683 
CHA0684 
CHAQ635 
CHA0686 
CHA0687 
CHA0688 
CHA0689 
CHA0690 
CHA0691 
CHA0692 
CHA0693 
CHA0694 
CHA0695 
CHA0696 
CHA0697 
CHA0698 
CHA0699 
CHA0700 
CHA0701 
CHA0702 
CHA07C3 
CHA0704 
CHA0705 
CHA0706 
CHAO 7  07 
CHA0708 
CHAU709 
CHA0710 
CNA0711 
CHA0712 
CHA0713 
CHA0714 
CHA0715 
CHA0716 
CHA0717 
CHA07  ]  8 
CHA0719 
CHA0720 
CHA0721 
CHA0722 
CHA0723 
CHA0724 
CHA0725 
CHA0726 
CHA0727 
CHA0728 
CHAQ729 
CHA0730 
CHA0731 
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CALL  BD0K1CL,X,P,DP.MIN,£> 

CALL  BDOKlCl,X,RO,i?RO,MIN,S) 

CALL  BDOICI.<L,X,Z,&Z  MM,4) 

PRINT  901,CNN,PRES$LNM),PULSE1CNN),PULSE2CNN), 

1  PULSE3CNN)/AMT0T,PULSE4CNN)/TENT0T,NN*1 >4) 

901  F0RMATC1X,2(V», I3,5D11.3,  '/•)/) 

IFCDABSCT-AC5)) .LT.l.D-6)  PRINT  911 , NERRP, ERRP 
911  F0RMAT(//1X,'NERRP,ERRP*',I5,D15.5/) 

RETURN 

7001  CONTINUE 

PRINT  7101,  I , ROAV , UAV , DRO ( I )  ,  DUC I ) ,  EC  I ) ,  EP , ZNEW , ZNEH-1 . DO , EPI 
7101  FORMAT!// IX, 'FROM  CYCEUL .  NEOATIVE  EP.  IN  CELL  I-»,I6// 

1  IX, 1 ROAV, UAV , DROC I ) , DUC I ) * • , 4D1S . 8// 

2  IX, ' ECI) , EP,ZNEW,ZN£W-1, EPI* 1 , 5D1 A .6//) 

CALL  PRINT 

1  (L,X,U,P,RO,G,E, DU, DP, DRO, DO, DXSI,MIN, 

2  U3, PS, UIDOT , PI DOT , 

*  FIMZ,ZMDOT , 

3  TENA, FIRO, FIM, FIE,GIP,VOL,Z,  DZ) 

CALL  SOFC 'CYCEUL  7001,  NEGATIVE  EP') 

RETURN 

- 5UlROUTTRF~S5TffE - SAfXeT 

1  (L,X,U,P,RO,G,E,DU, DP, DRO, DG, DXSI ,MIN, 

2  US, PS, UIDOT , PIDOT , 

X  FIMZ,ZMDOT , 

3  TENA, FIRO, FIM, FIE, GIP, VOL, Z,DZ) 

IMPLICIT  REALX8C A-H, Q-Z, ♦ ) 

DIMENSION  X(L),U(L),P(L),  RO(L),G(L),  E(L) ,  DU(l. ) ,  DP(  L)  ,  PROC  L ) , 


CHA0732 

CHA0733 

CHA0734 

CHA0735 

CHAQ736 

CHA0737 

CHA0738 

CHA0739 

CHA0740 

CHA0741 

CHA0742 

CHA0743 

CHA0744 

CHA0745 

CHAQ746 

CHA0747 

CHA0748 

CHA0749 

CHA0750 

CHA0751 

CHA0752 


CHA0755 
CHA0756 
CHA0757 
CHA0758 
CHA0759 
CHA0760 

1  DG( L ) , DXSI ( L ) , MINC  L ) ,  CHA0761 

2  USC  L ) , PSC  L ) , UIDOT ( L ) , PIDOTCL )  CHA0762 

3  ,TENA(L),FIRO(L),FIM(L),FIE(L)  CHA0763 

4  , GIP( L ) , VOL( L ) , Z( L) , DZ( L )  CHA0764 

5  , FIMZC  D , ZMDOTC  L )  CHA0765 

COMMON  /AB/AC50)  CHA0766 

EQUIVALENCE  (LL,A(2)),(T,A(3)),(DT,A!4)),( NCYC , A ( I 2 ) )  CHA0767 

EQUIVALENCE  (UGAL,AC15)J  CHA0768 

COMMON  /GAM/GAMA , NG, MU2, G1 , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , G1 0 , G1 1  CHA076 9 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23  CHA0770 

2  ,G24,G25.G26,G27,G28,G29,G30,G31,G32,G33,G34,G35  CHA0771 

REAL*8  NG,MU2  CHA0772 

COMMON/DETO/QDET , PCJDET , RCJDET , UCJDET , DCJDET , PODET , ROODET ,  CHA0773 

1  RATE, TEMPC  CHA0774 

C0MM0N/DIFFUS/U2, P2, R02, ARW  CHA0775 

CXXXXXXXKKXXXXi.XKKXXK*XK**X*XKXKXXXX*XXXXXXKXXX*KXKKXXXKXX*KXXKXXX*XXXXKCHA0776 
C  CHA0777 

C  RIGID  B.C.  AT  1=2  CHA0778 

C  CHA0779 


VC1)=-UC2.) 

PC1)=P(2) 

G(1)=G<2) 

R0(I)=R0C2) 

ZCD--ZC2) 

DUC 1)=DUC 2) 

DPC 1 ) =-DP( 2) 

DGC1)=-DGC2) 

DR0C1)*-PR0C2) 

DXSIC1 J=UXSIC2) 

C 

C  OUTFLOW  B.C.  AT  1=1 
C 

UCL)='JCin  +  DU(LL)/2.D0 
PCD *PC  ll  )+DPC  LD/2 .  DO 
ROC  L ) =ROC  LL )  +  DRO( LL  )/2 . DO 
GCL)=GCLL)+DGCLL)/2.D0 
ZCL)=ZCLD  +  DZCLD/2.D0 
DUC  L) =0 . 

DPC  L )  =  0 . 

DGC  L )  =  0 . 

DRQCL)=0. 

DZCL)  =  0  . 

DXSIC L )  =  DXSIC  LL  ) 


CHA0780 

CHA0781 

CHA0782 

CHAQ783 

CHA0784 

CHA0785 

CHA0786 

CHA0787 

CHA0788 

CHA0789 

CHA0790 

CHA0791 

CHA0792 

CHA0793 

CHA0794 

CHA0795 

CHA0796 

CHA0797 

CHA0798 

CHA0799 

CHA0800 

CHA0801 

CHA0802 

CHA0803 


PILE i  CHARGEP  FORTRAN  A1 


i-kp 
Ml) 


I'VfliMUlJL'J 


IMPLICIT  REALX8<A-H,0-Z,8)  *  CHA0808 

DIMENSION  X<L),V<L),DV<L),MIN<L)  CHA0809 

COMMON  /AB/AC50)  CHA0810 

EQUIVALENCE  ( LL , AC2) ) , (KEYMON, A<11 ) )  CHA08I1 

COMMON  /DRAW/GGDELX, GQDELY, UMIN, UMAX, PMIN, PKAX, ROMIN, ROMAX  CHA08I2 

1  ,XMIN,XMAX,SMIN,SMAX,IVERSA  CHA0813 

COMMON  /M0NIT/CASEAV<4),NC14<4),NF16<6),  CHA0814 

1  NMONU <  4  ) , NMONP <  4 ) , NMONRO <  4  ) »  NMONZ ( 4 )  CHA0815 

DIMENSION  NM0NV<4,4)  CHA0816 

EQUIVALENCE  ( NMONV (1,1), NMONU ( 1 ) )  CHA0817 

DIMENSION  NAMEV(A)  CHA08I8 

DATA  NAMEV/'U',  'P«,  'RO*,  *Z'/  CHA0819 

DATA  EPS/1. D-9/  CHA0820 

CXXXXXXXXXXXXXKXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXCHA0821 


•nm'M'jh 


GO  TO  <1,2,3, 4),  NV 
AMIDA-<UMAX-UMINJXX2 
GO  TO  9 

AMIDA*<PMAX-PMIN)XX2 
GO  TO  9 

AMIDA*(R0MAX-R0MIN)XX2 

GO  TO  9 

AMI DA *1. DO 

GO  TO  9 

CONTINUE 

AMIDA-AMIDAXEPSXX2 

E?*;A*DSQRT<AMIDA) 

DO  29  1*2, LL 
ICAT*0 

IF<DABS<DV<I )) .LE.EPSA)  DV<I)*0. 
IF<DV<I).EQ.O.)  GO  TO  29 
VLEFT =V( I)-0 . 5D0XDV< I ) 

VRIGHT =V<I)  +  0.5D0XDV<I) 

VM=V(I-1) 

VP=V(I+D 

SIGN=(VP-V(I) )X( V<I )-VM) 

IF< SIGN . GT . -AMI DA)  GO  TO  22 
DV< I ) =0 . 

ICAT  = 1 
GO  TO  20 
CONTINUE 

SIGN=( VP-VM)xDV(I) 

IF( SIGN . GT . -AMI DA )  GO  TO  24 
DV< I ) =0 . 5D0X( VP-VM) 

VLEFT =V<I)-0.5D0XDV<I) 

VRIGHT=V(I)+0.5D0*DV(I) 

ICAT=2 

SIGN=(VLEFT-VM)*DVC1) 

IF( SIGN . GT . -AMI DA)  GO  TO  26 
VLEFT=VM 

VRIGHT=2.D0XV(I) -VLEFT 
DVC I ) =VRIGHT-VL  EFT 
ICAT-3 

SIGN=(VP-VRIGHT)XDV(I) 

IF  (SIGN . GT . -AMIDA )  GO  TO  28 
VRIGHT  =  VP 

VLEFT  =  2 . DOXVC I ) -VRIGHT 

DV(I)=VRIGHT-VLEFT 

ICAT=3 

IF(DABS(DV(I)) .LE.0.5D0XDABSCVP-VM))  GO  TO  31 
DV(I)=0 .5D0X(VP-VM) 

ICAT=4 

CONTINUE 

CONTINUE 

IF  (DABSCDV(I)) .GT.EPSA)  GO  TO  40 
DV( I ) =0 . 

CONTINUE 

IF  (ICAT.GT.O)  NMONV ( ICAT , NV ) = NMONV (ICAT,NV)+I 
CONTINUE 
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CHA0831 

CHA0832 

CHA0833 
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CHA0836 

CHA0837 

CHA0838 

CHA0839 

CHA0840 

CHA0841 

CHA0842 

CHA0843 

CHA0844 

CHA0845 

CHA0846 

CHA0847 

CHA0848 

CHA0849 

CHA0850 

CHA0851 

CHA0852 

CHA0853 

CHA0854 

CHA0855 

CHA0856 

CHA0857 

CHA0858 

CHA0859 

CHA0860 

CHA0861 

CHA0862 

CHA0863 

CHA0864 

CHA0865 

CHA0866 

CHA0867 

CHA0868 

CHA0869 

CHA0870 

CHA0871 

CHA0872 

CHA0873 

CHA0874 

CHA0877 
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RETURN 


CHA0878 


3UBKUU I IRC  UUULCIL, A, 

IMPLICIT  REALX8C A-H, 0~Z, ♦ ) 
DIMENSION  XCL),VCL),DVCL),MINCL) 
COMMON  /AB/AC50) 

EQUIVALENCE  CLL,AC2)) 


J>COL£ 


CHA0880 

CHA088I 

CHA0882 

CHA0883 

CHA0884 


CXKXKKKXXMXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA0885 


DO  I  1*2, LL 

IM»I-1 

IP-I+1 

DV(I)*0.5D0X(V(IP)-V(IM)) 

CONTINUE 

RETURN 

END 


L  (L,X,U,P,RO,G,E,DU,DP, DRO, DG, DXSI ,MIN,  '  ^ 

1  US, PS, UI DOT , PIDOT , 

*  FIMZ, ZMDOT, 

S  TENA, FIRO, FIM, FI E, GIF, VOL, Z, DZ) 

IMPLICIT  REAL*8CA-H,0-Z> ♦  ) 

DIMENSION  XCL),UCL),PCL),ROCL),GCL),ECL),DUCL),DPCL), DRO(L ) , 
L  DOC  L ) , DXSI( L ) ,MIN( L ) , 

2  USC L ) , PSC L ) , UIDOTC L ) , PIDOTC L ) 

5  , TENACL ) , FIROC  L ) , FIMC  L ) , FIE( L ) 

»  , GIPCL) , VOLCL ) , Z( L ) , DZ( L ) 

5  , FIMZC  L ) , ZMDOT  CL) 

COMMON  /TOT/AMTOT, ETOT, EKTOT, EPTOT, TENTOT 

COMMON  /STEPO/UL , PL , ROL , OL , UR , PR , ROR, GR, USTAR, PSTAR, 

L  RSTARL , RSTARR, GSTARL , GSTARR, 

2  CL , CR, CSTARL , CSTARR, SL , SR , WL , WR, UWC  6 ) 


CHA0886 

CHA0887 

CHA0888 

CHA0889 

CHA0890 

CHA0891 

CHA0892 


CHAUBVO 

CHA0894 

CHA0895 

CHA0896 

CHA0897 

CHA0898 

CHA0899 

CHA0900 

CHA0901 

CHA0902 

CHA0903 

CHA0904 

CHA0905 

CHA0906 

CHA0907 

CHA0908 


, LAMDAL , LAMDAR, RATEL , RATER, TEMPL ,TEMPR,TEMPSL , TEMPSR  CHA0909 


4  , ZL , ZR, ZSTARL , ZSTARR, NFLUX, HEL  EML , HELEMR 

REALX8  LAMDAL, LAMDAR 
LOGICAL  HELEML, HELEMR 
COMMON  /AB/AC50) 

EQUIVALENCE  ( LL , AC2) ) . CT, AC3) ) , CNCYC, AC  12) ) , C DT, AC4 ) ) 
EQUIVALENCE  <UGAL,A(15)) 

COMMON/ D I FFUS/ U2, P2, R02, ARM 

COMMON/ DETO/QDET, PCJDET, RCJDET, UCJDET, DCJDET, PODET, ROODET, 
1  RATE, TEMPC 

COMMON  /CAM/GAMA, NG, MU2, G1 ,G2,G3, G4 ,G5,G6 , G7 ,G8, G9, GIO, Gil 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2  , G24 , G25 , G26 , G27 , G28 , G29 , G30 , G31 , G32 , G33 , G34 , G35 
REALX8  NG, MU2 

COMMON  /MONIT/CASEAVC  4 ) , NCI 4 C 4 ) , NF16 ( 6 ) , 

1  NMONUC  4 ) , NMONP C  4 ) , NMONRO  C  4 ) , NMONZC  4 ) 

DIMENSION  CA5AV1C4) 

LOGICAL  FULLPR 


CHA0910 

CHA0911 

CHA0912 

CHA0913 

CHA0914 

CHA0915 

CHAQ916 

CHA0917 

CHA0918 

CHA0919 

CHA0920 

CHA0921 

CHA0922 

CHA0923 

CHA0924 

CHA0925 

CHA0926 


CXXXXXXXXXXXXXXXXXMXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA0927 


FULLPR3 . TRUE . 

PRINT  1 
FORMAT ( 1H1 ) 

PRINT  2,  T , DT , NCYC 
FORMAT(1X,10X, 'RESULTS  AT  T= ' , Dll . 5, 5X, ' DT= • , D1 1 . 5, 5X, ' NCYC= ' 
1  15//) 

PRINT  3,  AMTOT, ETOT, EKTOT, EPTOT, TENTOT 

FORMAT ( IX, • AMTOT= ' , D20 . 14 , 2X , ' ETOT, EKTOT, EPTOT= ' , 3D22 .14/ 

1  IX, *TENT0T=',D21.14//) 

FORMAT  MX.*  T»  •  Y  •  »  II  ».*  P 


FORMAT ( IX, 


FORMAT (IX, 


U 

G 

DP 

DZ') 

US 

FIMZ 

TEMP 

ENTRO 


L  '  ZMDOT  FIMZ  ' 

l  '  AMDOTN  TEMP 

5  '  AMACH  ENTRO  ') 

FORMAT (IX) 

IF  (UGAL.NE.O.)  PRINT  6,  UGAL 

FORMATC/llX, 'INITIAL  VELOCITY  CORRESPONDS  TO 

DO  10  1=1, L 

IF  (MODC I , 10) . NE . 1 )  GO  TO  11 


1  PS 
1  AMDOT 
ENTALP 


UGAL= ' , D15 . 6/) 


CHA0928 
CHA0929 
CHA0930 
CHA0931 
CHAO 9 32 
CHA0933 
CHA0934 
CHA0935 
CHA0936 
CHA0937 
CHA0938 
CHA0939 
CHA0940 
CHA0941 
CHA0942 
CHAQ943 
CHA0944 
CHA0945 
CHA0946 
CHA0947 
CHA0948 
CHA0949 


FILE*  CHARGE?  FORTRAN  A1 


PRINT  3 
PRINT  4 
PRINT  44 
PRINT  5 

11  CONTINUE 

PRINT  12,I,X<I),U<I),P(I),R0<I>,G<I),Z(I),DUU),DPa),DR0<I>, 

1  BG( I ) ,  DZ( I ) 

12  F0RMAT(1X,I3,6D12.5,5D11 .4) 

ENTRO*P( I )/RO( I )X*GAMA 
IF(.NOT.FULLPR)  GO  TO  131 
IF(I.EQ.l)  GO  TO  131 
IM“I-1 

UL*UC IH)+0 . 5*DU( IM) 

PL«PCIM)+0.5*DP(IM) 

ROL  *RO( IM)+0 . 5*DRO( IM) 

GL«G(IM)+0.5*DG(IM) 

Cl-Gl/ROL 

ZL*Z( IM)+0 . 5*DZ( IM) 

IFCZL.LT.O.)  ZL*0 . 

UR»U<I)-0.5*DU(I) 

PR«P( I)-Q . 5*DPC I ) 

GR«G(I)-0.5*DG(I) 

ROR=RO(I)-0.5XDRO(I) 

CR»GR/ROR 

ZR*Z( I ) -0 . 5*DZ( I ) 

IFCZR.LT.O.)  ZR=0 . 

IFCPL . LE. 0 . )  PLS1 . D-8 
IFCPR.LE.Q.)  PR*1 . D-8 
CALL  RIEMANCL,I,MIN) 

XI»X(I) 

RSTAR-RSTARL 

IF<USTAR.LT.O. )  RSTAR*RSTARR 
ZSTAR-ZL 

IFtUSTAR .  LT . 0 . )  ZSTAR=ZR 
AMACH=U5TAR/DSQRT(GAMAXPSTAR/RSTAR) 

AMDOT»RSTARXU$TARXCROSS(XI) 

IFCI.NE.2)  GO  TO  132 
AMDOTO*AMDQT 

I F( DABS (AMDOTO ) . LT . 1 . D-12)  AMDOTO=1.DO 
132  CONTINUE 

AMDOTN=AMDGT/ AMDOTO 

ENTALP=(GAMA/(GAMA-1.D0))*PSTAR/RSTAR+Q.5D0XUSTAR**2+QDET*ZSTAR 
ARW=1 . DO 

TEMP=PSTAR/  (  RSTARXARW) 

PRINT  13,US(I),PS(I), 

1  ZMDOTCI ) , FIMZC I ) , AMDOT, AMDOTN, TEMP , ENTALP, AMACH , ENTRO 

13  FORMAT  (4X,12X,5D12.5,6D11.4) 

131  CONTINUE 

10  CONTINUE 
C  JOB  STATISTICS 
DO  40  1=1,4 
CASAV1(I)=0 . 

IF  (NC14(I) .NE.O)  CASAV1 ( I )=CA5EAV( I )/DFL0AT(NC14( I ) ) 

40  CONTINUE 
PRINT  30 

30  F0RMAT(///1X,10('X«),3X, 'JOB  STATISTICS  1 , 3X, 1 OC • *' )// ) 

PRINT  31, (NC 14(1), 1=1, 4) 

31  FORMAT (IX, 1  NO .  QF  VARIOUS  CASES  IN  RIEMAN  SOLVER  NC14(NCASE)= ' , 

1  4110) 

PRINT  301,  (CASAVl(I) ,1=1,4) 

301  F0RMAT(/1X, 'AVERAGE  NUMBER  OF  ITERATIONS  IN  RIEMAN  SOLVER', 

1  IX,'  CAS AVI ( NCASE) = ',4(F6 .2,4X)) 

PRINT  32,(NF16(I),I=1,6) 

32  F0RMAT(/1X, 'NO.  OF  VARIOUS  FLUX  CASES  NF16(NFLUX)=' ,6110) 
ICATQ=4 

PRINT  33, (NMONU( I ) , 1=1 , ICATO) , (NMONP(I ) , 1=1, ICATO ) , 

1  (NMONRO(I), 1=1, ICATO) , (NMONZ(I), 1=1, ICATO) 

33  F0RMAT(/1X, 'NO.  OF  MONOTONICITY  INTERVENTIONS  FOR  EACH  VAR.', 

1  IX, 'IN  EACH  CATEGORY.'/ 

1  IX, 'NMONU  ( ICAT ) = ' , 4110/ 

1  IX, 'NMONP  ( ICAT ) =  ' , 411 0/ 

1  IX, 'NMONRO(ICAT)  =  ',4HO/ 
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CHA1008 
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CHA1012 

CHA1013 

CHA1014 

CHA1015 
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CHA1019 
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CHA1023 

CHA1024 

CHA1025 

CHA1026 

CHA1027 

CHA1028 
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RETURN 


IX# 'NHONZ  (ICAT)*' #4110/) 


SUBROUTINE  RIEMANCL, I, MIN)  0\ 

IMPLICIT  REALX8CA-H, O-Z, $)  n ' 

DIMENSION  MIN(L) 

COMMON  /STEPO/UL ,  Pi.  #  ROL , GL ,  UR ,  PR ,  ROR ,  GR ,  USTAR,  PSTAR, 

1  RSTARL  #  RSTARR#GSTARL , GSTARR, 

2  CL , CR>  CSTARL #  CSTARR#  SL  , SR, ML ,  HR, UW(6 ) 


SOF 


tflEM  AN' 


CHA1030 

CHA1031 


IMPLICIT  REAL*8(A-H»  O-Z# ♦)  *  n  CHA1035 

DIMENSION  ISTOP(l)  CHA1036 

PRINT  1#  ISTOP  CHAi.037 

1  FORMAT (//IX, 5HXXX, 2X, 20A4, 3X,  3HXXX///)  CHA1038 

PRINT  1  CNA1Q39 

XX—l.DO  CHA1D40 

YY*DSQRT(XX)  CHA1041 

STOP  CHA1042 

_ EHB _ r - QULinAS- 

SUBROUTINE  RIEMANCL,  I, MIN)  R\F  MAN  CHA1310 

IMPLICIT  REALX8(A-H,0-Z,$)  n n  CHA1311 

DIMENSION  MIN(L)  CHA1312 

COMMON  /STEPO/UL, PL, ROL, GL, UR, PR, ROR, GR, USTAR, PSTAR,  CHA1313 

1  RSTARL, RSTARR,GSTARL, GSTARR,  CHA1314 

2  CL, CR, CSTARL, CSTARR, SL, SR, WL, HR, UW<6)  CHA1315 

3  #  LAMDAL , LAMDAR, RAT EL , RATER, TEMPL , TEMPR, TEMPSL , TEMPSR  CHA1316 

4  , ZL , ZR, ZST ARL , ZSTARR , NFLUX#  HEL  EML , HELEMR  CHAI3I7 

REALK8  LAMDAL, LAMDAR  CHA1318 

LOGICAL  HELEML, HELEMR  CHA13I9 

COMMON  /STEPI/DUIDT , DPIDT , DGIDTL , DGIDTR, DRIDTL , DRIDTR  CHAI320 

2  , ASTARL, ASTARR, LAMDSL , LAMDSR, DSDAL#  DSDAR, DZDAL, DZDAR  CHA1321 

3  , RAT  #  SH  CHA1322 

4  , BETACL , BETACR, DSDASL , DSDASR, DZDASL , DZDASR  CHA1323 

REALX8  LAMDSL, LAMDSR, DSDAL, DSDAR, DZDAL, DZDAR  CHA1324 

COMMON  /DRAW/GODELX,GODELT, UMIN, UMAX, PMIN, PMAX, ROMIN, ROMAX  CHA1325 

l  , XMIN,XMAX, SMIN, SMAX, I VERSA  CHA1326 

COMMON  /GAM/GAMA , NG, MU2,G1 , G2, G3, G4, G5,G6 ,G7 , G8,G9,G10, Gil  CHA1 327 

1  ,G12#G13,G14,G15,G16,GI7,G18,G19,GZ0,G21,G22,G23  CHA1328 

2  , G24 , G25, G26 , G27 , G28 , G29 , G30 , G31 , G32 , G35, G34, G35  CHA1329 

REAL *8  NG, MU2  CHA1330 

COMMON  /AB/AC50)  CHA1331 

COMMON  /M0NIT/CASEAVC4) , NC14(4),NF16(6),  CHA1332 

1  NM0NU(4 ) , NM0NPC4? , NM0NR0(4) , NMONZC  4)  CHA1333 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA1334 
DATA  NMAX/63Z  CHA1335 

DATA  EPS/1. D-8/  CHA1336 

DATA  NTRY/O/  CHA1337 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA1338 
UW(6)=1.D20  CHA1339 

WL  =  0.  CHA1340 

WR-0.  CHA1341 

ZETAL=PLXXG8  CHA1342 

ZETAR=PRXXG8  CHA1343 

CLG=CL/GAMA  CHA1344 

CRG=CR/GAMA  CHA1345 

ZSTARL  =  ZL  CHA1346 

ZSTARR=ZR  CHA1347 

IF  (ZETAL .LT.ZETAR)  GO  TO  102  CHA1348 

C  LEFT  PRESSURE  IS  HIGHER  CHA1349 

101  CONTINUE  CHA1350 

EVERR=(PL-PRVPR  CHA1351 

USR=UR+CRGXEVERR/DSQRT(1 . DQ+G6XEVERR)  CHA1352 

SRR=USR  CHA1353 

UEL=UL-G7XCLX(ZETAR-ZETAL)/ZETAL  CHA1354 

SLL=UEL  CHA1355 

NL  =2  CHA1356 

NR=2  CHA1357 

IF  CUSR.GE.UU  NL  =  1  CHA1358 

IF  (UEI..LE.UR)  NR=1  CHA1359 

IF  ( DABSC  EVERR) . LT . EPS )  GO  TO  100  CHA1360 

IF  (NL.EQ.2.AND.NR.EQ.1)  GO  TO  7001  CHA1361 

GO  TO  100  CHA13A2 

C  RIGHT  PRESSURE  IS  HIGHER  CHA1363 

102  CONTINUE  CHA1364 

EVERL=(PR-PL)/PL  CHA1365 

USL=UL-CLGXEVERL/DSQRT(1 .D0+G6XEVERL)  CHA1366 

SLL=USL  CHA1367 


GO  TO  100 

C  RIGHT  PRESSURE  IS  HIGHER 
102  CONTINUE 

EVERL=(PR-PL)/PL 

USL=UL-CLGXEVERL/DSQRT(1 .D0+G6XEVERL) 
SLL=USL 


UER=UR+G'/  XCRX (  ZETAL  -ZETAR)/ZETAR 


CHA1338 


FILE*  CHARGEP 
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SRR=UER 

NL»2 

NR*2 

IF  (UER.GE.UD  NLaI 
IF  CUSL.LS.UR)  NR=1 
IF  (DABS(EVERL) .LT.EPS)  GO  TO  100 
IF  (NL.EQ.l . AND . NR . EQ . 2)  GO  TO  7001 
GO  TO  100 
100  CONTINUE 

IF  (Nl.EQ.l.AND.NR.EQ.2)  NCASE=1 
IF  (NL . EQ . 2 . AND . NR . EQ . 2)  NCASE=2 
IF  (NL . EQ . 2 .AND . NR . EQ . 1 )  NCASE=3 
IF  (NL.EQ.l.AND.NR.EQ.l)  NCASE=4 

IF(DABS(PL-PR)+DAbS(UL-UR) . LT . EPS#( PMAX-UMIN) )  NCASE=4 
UMIDA=EPS*DMAX1(CL,CR) 

DUDZL*~G7*CL/ZETAl 
DUDZR3  G7XCR/ZETAR 

ZETA=(-(UR-UL)+Z£TARKDUDZR-ZETAL»DUDZL)/(DUDZR-DUDZL) 

IF  (ZETA.LE.O.)  GO  TO  7002 

N-0 

GO  TO  (1,2, 3, 4),  NCASE 
C  THE  CASE  ES 

1  ITYPE=NCASE 
HELEML3. FALSE. 

HELEMR3 .TRUE. 

11  N=N+1 

IF  (N.GT.NMAX)  GO  TO  7003 
ZETAF=ZETA 

UEL3UL-G7*CL*(Z£TAF-ZETAL)/ZETAL 

PPR3(ZETAF/ZETAR)**NG 

EVERR=PPR-i . DO 

SQRR=DSQRT(1 .D0*G6*EVERR1 

USR=UR+CRG*EVERR/SQRR 

DU=UEL-USR 

IF  (DABS (DU) . LE.UMIDA)  GO  TO  10 

DUDZR  =  NG*CRGX(r'PR/ZETAF)*(l .  D0+G9*EVERR  VSQRR*«3 
ZETA=ZETAF+DU/ ( DUDZR-DUDZL ) 

GO  TO  11 
10  CONTINUE 

USTAR=(UEL+USR)/2.D0 
IF(DABSCUSTAR) .LT.EPSXUMAX)  USTAR  =  0  . 

PSTAR=PPR*PR 

CSTARL=CL+(UL-USTAR)/G7 

RSTARl=OAMAXPSTAR/CSTARL**2 

^TARI  *R^TAPI 

C  EQU.  NO.  69.01  OF  THE  BOOK  BY  COURANT-FRIEDRICHS . 
WWR=G11*(USTAR-UR)XR0R 
WR=WWR+DSQRT(GRX*2+WWRX*2) 
RSTARR=RORXWR/(WR-ROR*(USTAR-UR)) 
GSTARR=DSQRT(GAMAXPSTARXRSTARR) 

CSTARR=GSTARR/RSTARR 

WRE=WR/ROR+UR 

UW(1)=UL-CL 

UW( 2) =USTAR-CSTARL 

iJW(  3 )  =USTAR 

UW(  4 ) =WRE 

UW( 5 ) =WRE 

GO  TO  5 

C  THE  CASE  SS 

2  ITYPE=NCASE 
HEL.  EML  =  .  TRUE  . 

HELEMR= . TRUE . 

21  N=N+1 

IF  (N.GT.NMAX)  GO  TO  7003 

ZETAF=ZETA 

PF=ZETAF**NG 


<y.  a  r,f' 


PPL=PF/PL 

PPR=PF/PR 

EVERL  =PPL-1 . DO 

EVERR=PPR-1 . DO 

SQRL "D3QRT ( 1 . D0+G6*EVERL  ) 

SQRR=DSQRT ( 1 . D0+G6XEVERR) 
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CHA1330 

CHA1381 

CHA1382 

CHA1383 

CHA1384 

CHA1385 

CHA1386 

CHA13S7 

CHA1388 

CHA1389 

CHA1390 

CHA1391 

CHA1392 

CHA1393 

CHA1394 

CHA1395 

CHA1396 

CHA1397 

CHA1398 

CHA1399 

CHA1400 

CHA1401 

CHA1402 

CHA1403 

CHA1404 

C.HA1405 

CHA1406 

CHA1407 

CHA1408 

CHA1409 

CH«1410 

CHA1411 

CHA1412 

CHA1413 

CHA1414 

CHA14J.5 

CHA1416 

CHA1417 

CHA1418 

CHA1419 

CHA1420 

CHA1421 

CHA1422 

CHA1423 

CHA1424 

CHA1425 

CHA1426 

CHA1427 

CHA1428 

CHA1429 

CHA1430 

CHA1431 

CHA1432 

CHA1433 

CHA1434 

CHA1435 

CHA1436 

CHA1437 

CHA1438 

CHA1439 

CHA1440 


m £h£i ^ 


CHARGER  FORTRAN  A1 


USL"UL-*CL0*EVERL/StJRL 

USR«UR+CRG*EVERR/SQRR 

DU-USL-USR 

IF  (DASS(DU) .LE.UMIDA)  GO  TO  20 
DUDZL—NO*CLO*(PPL/ZETAF)X(1.DO+G9*EVERL)/SQRLXX3 
DUDZR*  NG*CRGX(PPR/ZETAF>*C1.DO+G9X£VERR)/SORRXX3 
ZETA-ZETAF+DU/ <  DUDZR-DUDZL ) 

GO  TO  21 
20  CONTINUE 

USTAR*(USL+USR)/'2 .  DO 

IF( DABStUSTAP.)  ,  IT . EPSXUMAX)  USTAR=0. 

PSTAR«(PPLKPL+PPRXPR)/2.D0 

NNft«GllX(USTAR-UR)XROR 

HR«WWR+DSQRT  (  GRxx2+WHRX*2  ) 

WWL*-G11X(USTAR-UL)XR0L 

WL*HWL+D$QRT(GLXX2+WWL*X2) 

RSTARL»ROLXWL/< WL+ROLXC  USTAR-UL ) ) 
RSTARR«RORXWR/(HR-RORX(USTAR-UR)) 
G$TARL*DSQRT<GAMAXPSTAR*RSTARL ) 
GSTARR*DSQRT<GAMAXPSTAR*RSTARR) 

CSTARL*GSTARL/RS  TARL 

CST  ARR*GST  ARR/RST  ARR 

HLE*-WL/RQL+UL 

WRE“WR/ROR+UR 

UW(1)*WLE 

UW(2)»WLE 

UW(3)»USTAR 

UH(4)*WRE 

UH(5)*WRE 

GO  TO  5 

C  THE  CASE  SE 

3  ITYf’E*NCASE 
HELEML1*  .TRUE. 

HELEMR=. FALSE. 

31  N=N+1 

IF  (N.GT.NMAX)  GO  TO  7003 
ZETAF3ZETA 

UER=UR+G7*CRX(ZETAF-ZETAR)/ZETAR 

PPL  =  (ZETAF/ZETAUXXNG 

EVERL=PPL~1.D0 

SQRL  =  DSQRT(1  .D0+G6XEVERU 

USL3UL-CLGXEVERL/SQRL 

DU-USL-UER 

IF  (DAPSCDU) .LE.UMIDA)  GO  TO  30 

DUDZL=-NGXCIGX(PPI/'ZETAF)*<1.D0+G9*EVERL)/5QRLX*3 
ZETA=ZETAF+DU/ ( DUDZR-DUDZL ) 

GO  TO  31 
30  CONTINUE 

USTAR-(USL+UER)/2.DG 
IF(DABS(USTAR) .LT.EPSXUMAX)  USTAR=0. 

PSTAR=PPLXPL 

CSTARR=CR-( UR-USTAR )/G7 

RSTARR=GAMAXPSTAR/CSTARRX*2 

GSTARR=CSTARR*RSTARR 

WWL=-G11X(USTAR-UL)XR0L 

WL-WWL+DSQRT(GLXX2+WWLX*2) 

WLE’-WL/ROL  +  'JL 

RSTARL=R0L*NL/'WL+R0LX(USTAR-UL)1 

GSTARL=DSQRT(GAMAXPSTARXRSTARL ) 

CSTARL=GSTARL/RSTARL 

UW(1)=WLE 

UWC2)=WLE 

UW(3)=USTAR 

UW(4)=USTAR+CSTARR 

UW( 5) =UR+CR 

GO  TO  5 

C  THE  CASE  EE 

4  ITYPEsNCASE 
HELEML-. FALSE. 

HELEMR=. FALSE. 

PSTAR=ZETAXXNG 

USTAR=UL-G7XCLX(ZETA-ZETAU/ZETAL 
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CHA1446 
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CHA1454 

CKA1455 

CHA1456 

CHA1457 

CHA1458 

CHA1459 

CHA146Q 

CHA1461 

CHA1462 

CHA1463 

CHA1464 

CHA1465 

CHA1466 

CHA1467 

CHA1468 

CHA1469 

CHA1470 

CHA1471 

CHA1472 

CHA1473 

CHA1474 

CHA1475 

CHA1476 

CHA1477 

CHA1478 

CHA1479 

CHA1480 

CHA1481 

CHA1482 

CHA1483 

CHAI484 

CHA1485 

CHA1486 

CHA1487 

CHA1488 

CHA1489 

CHAI490 

CHA1491 

CHA1492 

CHA1493 

CHA1494 

CHA1495 

CHA1496 

CHA1497 

CHA1498 

CHA1499 

CHA1500 

CHA1501 

CHA1502 

CHA1503 

CHA1504 

CHA1505 

CrtAl  506 

CHA1507 

CHA1508 

CHA15Q9 

CHA1510 

CHA1511 

CHA1512 


FILE t  CHARGEP  FORTRAN  A1 


IF(DABSCUSTAR) . LT . EPSXUMAX)  USTAR«0.  CHA1513 

CSTARL*CL+(UL -USTAR )/G7  CHA1514 

CSTARR<*CR-(UR-USTAR)/G7  CHA1515 

R$TARL*GAMA*PSTAR/C$TARLX*2  CHA1516 

R$TARR*GAMAXPSTAR/C$TARRXX2  CHA1517 

GSTARL»RSTARLXCSTARL  CHA1518 

GSTARR«RSTARRXCSTARR  CHA1519 

UW(1)*UL~CL  CHA1520 

UW(2)°USTAR-CSTARL  CHA1521 

UW(5)*USTAR  CHA1522 

UW(4) “US7AR+CSTARR  CHA1523 

UN(5)«UR+CR  CHA1524 

N«1  CHA1525 

GO  TO  5  CHA1526 

5  CONTINUE  CHA1527 

DO  6  K*l,6  CHA1528 

NFLUX*K  CHA1529 

IF  (UW(K).GE.O.)  GO  TO  61  CHA1530 

6  CONTINUE  CHA1531 

NFLUX*6  CHA1532 

61  CONTINUE  CHA1533 

NC14(NCASE)«NC14CNCASE)+1  CHA1537 

CASEAV(NCASE)*CASEAV<NCASE)+DFLOATCN)  CHA1538 

NF16(NFLUX)aNF16(NFLUX)+l  CHA1539 

IF(NTRY.GE.2>G0  TO  666  CHA1540 

IF( I . NE. 2 . AND. I . NE . L )  GO  TO  666  CHA1541 

PRINT  667,1, NFLUX, NCASE, PL , UL , ROl , PR, UR, ROR, USTAR, PSTAR, RSTARL ,  CHA1542 

1  RSTARR, ( KK, UW<  KK) , KX=1 , 6 )  CHA1543 

66?  F0RMAT</1X, • I, NFLUX,NCASE= • , 3I5/1X, ' PL , UL , ROL , PR, UR, ROR* ' ,6D12 . 4/  CKA1544 

1  IX, 'USTAR, PSTAR, RSTARL, RSTARR= • , 4D13 . 4/  CHA1545 

2  IX, 'KK, UM(KK) a',6(I4,2X/D13.4)/)  CHA1546 

NTRYaNTRY+l  CHA1547 

666  CONTINUE  CHA1548 

RETURN  CHA1549 

7001  CONTINUE  CHA1550 

PRINT  7101,  PL , UL, PR, UR, ZETAL , ZETAR >  SLL, SRR, NL , NR, I  CHA1551 

7101  F0RMATC//1X, 'FROM  RIEMAN.  AN  IMPOSSIBLE  CASE  OF  EXPANSION/SHOCK '  CHA1552 

1  //IX, 8Pl,Ul,PR,UR=’,4D25.14//  CHA1553 

2  IX, ’ZETAL, ZETAR, SLL, SRRa’,4D25. 14//  CHA1554 

3  IX, 'NL, NR. I*', 3110//)  CHA1555 

CALL  SOFC7001*)  CHA1556 

7002  CONTINUE  CHA1557 

PRINT  7102,  ZETA , DUDZL ,  DUDZR, ZETAL , ZETAR ,PL,UL,PR,UR,N,  NCASE,  I  CHA1558 

7102  F0RMATI//1X, ’FROM  RIEMAN.  NEGATIVE  PRESSURE  AT  THE  INTERSECTION », CHA1559 

1  IX, ’OF  L  AND  R  EXPANSION  BRANCHES’//  CHA1560 

2  IX, ’IT  MEANS  THAT  A  CAVITATION  TENDS  TO  FORM.  THIS’,  CHA1561 

3  IX, ’POSSIBILITY  IS  EXCLUDED  IN  PRESENT  VERSION’//  CHA1562 

4  IX, ’ ZETA, DUDZL , DUCZR: ZETAL , ZETAR, PL , UL , PR, UR= ’ , 9D10 . 3//  CHA1563 

5  IX, *N, NCASE, 1=’, 3110//)  CHA1564 

CALL  SOF( ’ 7002 ’ )  CHA1565 

7003  CONTINUE  CHA1566 

PRINT  7103,  I , N, NCASE, DU, UMIDA, EPS, PL, UL, PR, UR,  CHA1567 

1  ZETA, ZETAF, ZETAL, ZETAR, DUDZL,DUDZR  CHA1568 

7103  FORMATC//1X, ’FROM  RIEMAN.  NUMBER  OF  ITERATIONS  EXCEEDED.’//  CHA1569 

1  IX, ’I, N, NCASE, DU,UMIDA,EPS=’,3I6,3D18.6//  CHA1570 

2  IX, ’PL,UL,PR.UR,ZETA,ZETAF=’,6D18.10//  CHA1571 

3  IX, ’ZETAL, ZETAR, DUDZL , DUDZR= » , 4D18 . 10//)  CHA1572 

CALL  SOFC7003’)  CHA1573 

RETURN  CHA1574 

END  CHA1575 


SUBROUTINE  MAGA  (  L ,  I ,  MIN )  Px 

IMPLICIT  REAL*8(A-H,0-Z,$)  1  ,no  ' 

DIMENSION  MINI L ) 

COMMON  /GAM/GAMA, NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10, Gil 

1  ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2  , G24 , G25, G26 , G27 , G28 , G29 , G30 , G31 , G32 , G33, G34 , G35 
REALXg  NG.MU2 

COMMON/ DETO/QDET , PCJDET , RCJDET , UCJDET, DCJDET, PODET, ROODET, 
1  RATE, TEMPC 

COMMON  /STEPO/UL, PL, ROL, GL, UR, PR, ROR, GR, USTAR, PSTAR, 

1  RSTARL, RSTARR, GSTARL , GSTARR, 
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CHA1578 

CHA1579 

CHA1580 

CHA1581 

CHA1582 

CHA1583 

CHA1584 

CHA1585 

CHA1586 

CHA1587 


CHARGEP  FORTRAN  A1 


PAGE 


Cl , CR, CSTARl , CSTARR, 5L , SR, WL , WR, UW! 6  >  CHA1588 

, L AMD A l ,  LAMDAR, RATEL , RATER, TFMPL ,  TEMPR, TEMPSL , TEMP SR  CHA1589 

CHA1590 


2 

3 

4  , ZL , ZR , ZSTARL , ZSTARR, NFLUX, HEL EML , HELEMR 
REALMS  LAMDAL , LAMDAR 
LOGICAL  HELEML, HELEMR 

COMMON  /STEP1/ DUIDT , DPIDT, DGIDTL , DGIDTR, DRIDTL , DRIDTR 

2  , ASTARL , ASTARR, LAMDSL , LAMDSR, DSDAL , DSDAR, DZDAL , DZDAR 

3  , RAT , SH 

4  , BETACl , BETACR , DSDASL , DSDASR , DZDASL , DZDASR 
REALMS  LAMDSL , LAMDSR, DSDAL , DSDAR, DZDAL , DZDAR 
COMMON  /GRADS/ DUDXIL , DPOXIL , DGDXIL , DRDXIL , DZDXIL , DSDXIL , 

1  DUDXIR, DPDXIR, DGDXIR, DRDXIR, DZDXIR, DSDXIR 

COMMON  /AB/A<50) 

REALMS  LU, LP, LRO, LLAMDA 
DATA  EPS/1. D-6/ 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA1603 
C  WE  HERE  SOLVE  FOR  THE  TIME-DERIVATIVES  ALONG  THE  CONTACT  SURFACE,  CHA1604 
NAMELY  DUIDT, DPIDT.  FROM  THESE  WE  ALSO  OBTAIN  THE  OTHER  CHA1635 

TIME-DERIVATIVES  (SEE  COMMON  /STEP1/).  CHA16G6 

WE  COMPUTE  THE  COEFFICIENTS  FOR  TWO  EQUATIONS  FOR  DUIDT, DPIDT.  THESECHA1607 
ARE  AALMDUIDT+BBLXDPIDT=DDL  CHA1608 

AARXDUIDT+B8RXDPIDT=DDR  CHA1609 


CHAI591 

CHA1592 

CHA1593 

CHA1594 

CHA1595 

CHA1596 

CHA1597 

CHA1598 

CKA1559 

CHA1600 

CHA1601 

CHA1602 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX*XXXXXXXXXXXXXKXCHA16) 0 


c 

c 

c 


IF(SH.LE.EPS)RAT=0. 


LEFT  SIDE  OF  CONTACT 


IF  ( .NOT.HELEML)  GO  TO  12 
11  CONTINUE 
C  LEFT  SHOCK 

DP=PSTAR-PL 

DU=USTAR-UL 

Z2=0 . 5D0/»'  PSTAR+MU2XPL  ) 

LU=DUX(0.5D0XR0L+MU2XZ2XGLXX2)-GLXX2/WL-WL 
LR0=-0 . 5DOXDP/ROL 
LP=-2.D0-MU2XZ2XDP 
AAL  =  2 . D0-Z2XCP 

BBL=Z2XDU+WL/GSTARLXX2+1 .DO/WL 
DDmUXDUDXIL+LROXDRDXIL  +  LPXOPnXIl 
DDL  =  DDL-WL*US‘f  ARxRAT/RSTARL 
1  +ULxRATX(-GAMAXPL/WL+DUXCGAMAXPLXMU2xZ2+0 . 5D0 ) ) 
GO  TO  10 
CONTINUE 
LEFT  RAREFACTION 

A1=DUDXIL+DPDXIL/GL 
BETA=GSTARL/Gl 
SQB=DSQRT (BETA) 

ASTARL=AI-(CI/(G15XSL))XDSDXILX(BETAXXG5-1 .DO) 
AAL=1 .DO 
BBL=1 .DO/GSTARL 
DDL—GSTARLXASTARL/SQB 
DSDAL=DSDXIL 


12 


DZDAL  =DZDXIL 
DSDASL~DSDXILXSQB 
DZDASL=DZDXILXSQB 
GEOM=RATXC (GAMA-1 . 


10 


C 

C 

C 


D0)*UL+2.DQ*CL)x 
1  (BETAXXG13-1 .DO)/ (ROLX( GAMA-3 . DO) ) 

1  -4.D0XRATXCIXCBETAXXG14-1 . DO)/(ROLX<3, 
ASTARL = AST ARL-GEOM 
EVER1 3  GSTARLXGEOM/SQB 
EVER2=-RAT*USTAR*CSTARL 
DDL  =DDl +EVER1  +  EVER2 
GO  TO  10 
CONTINUE 


DOXGAMA-5 . DO ) ) 


RIGHT  SIDE  OF  CONTACT 


IF  (  .NOT. HELEMR) 
21  CONTINUE 
C  RIGHT  SHOCK 

DP=PSTAR-PR 

DU=USTAR-UR 


GO  TO  22 
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Z2*0 .5D0/( PSTAR+MU2*PR ) 

LU*DU*( 0  5DQXR0R  MU2«Z2*0RXX2)+GRXX2/WR+HR 
LRO><  0.500MDP/R0R 
IP*-2.D0-MU2*Z2*DP 
AAM2.00-Z2MDP 

BBR»Z2*DU-WR/G5Y ARRXM2-1 . DO/NR 

dor*lu*dudxir+lroxdrdxir+lp*dpdxir 

DDR*DDR+WRxUSTAR«RAT/R$TARR 
1  ♦URXRATXt  GAMA*PR/WR+DU*<  0AMA*PR*MU2*Z2+O .500)) 


1 

2 


2 

3 

4 


CHA166Q 

CHA1661 

CHA1662 

CHA1663 

CHA1664 

CHA1665 

CHA1666 

CHA1667 

CHA1668 


00  TO  20 

CHA1669 

22  CONTINUE 

CHA1670 

C  RIGHT  RAREFACTION 

CHA1671 

Al-DUDXIR-DPDXIR/GR 

CHA1672 

BETA-OSTARR/GR 

CHA1673 

SQB*D5QRT( BETA) 

CHA1674 

ASTARR»A1-KCR/<G15XSR))*DSDXIRX(BETAXXG5-1.D0) 

CHA1675 

AAR*1 . DO 

CHA1676 

BBR*-1 . DO/OSTARR 

CHA1677 

ODR>OSTARRXASTARR/SQB 

CHA1678 

DSDAR” DSDXIR 

CHA1679 

DZDAR* DZDXIR 

CHA168C 

DSDASR-DSDXIRXSQB 

CHA1631 

DZDASR*DZDXIRXSQB 

CHA1682 

QEOM*RATW( - <  GAMA-1 . DO ) KUR+2 . D0XCR)X( BETAXKG1 3-1 . DO ) 

CHA1683 

1  /( R0RX(GAMA"3 . DO) ) 

CHA1684 

2  -4 . DOXRAT«CR*<  BETA«*G14-1 . DO )/ ( ROR*( 3 . DQ*GAMA-5 . DO ) ) 

CHA1685 

A5TARR“ASTARR+GE0M 

CHA1686 

EVER1*GSTARR*3£0M/SQB 

CHA1687 

EVER2-RATXUSTARXCSTARR 

CHA1688 

DDR-DDR+EVER1+EVER2 

CHA1689 

00  TO  20 

CHA1690 

20  CONTINUE 

CHA1691 

DET-AALXBBR-AARXBBL 

CHA1692 

DUIDT-(DDLKBBR'DDRXBBL)/DET 

CHA1693 

DPIDT*-(DDLXAAR-DDRXAAL)/DET 

CHA1694 

DRIDTL*DPIDT/CSTARLXXZ 

CHA1695 

DRIDTR*DPIDT/CSTARRXX2 

CHA1696 

RETURN 

CHA1697 

_ £N£ _ 

CHA1698 

— — in  iii  I'niii  imyKHir— —  i«n  — 

crraror? 

IMPLICIT  REALX8(A-H,0-Z,$)  r  L.  ^  *  C.  CHAI'OO 

DIMENSION  MINI L )  CHAI701 

COMMON  /AB/AI50)  CHA1702 

EQUIVALENCE  (DT, AC4) ) , (NCYC, A( 12) )  CHA1793 

COMMON  / GAM/ GAMA, NG, MU2, G1 , G?  G3 , G4 , G5, G6 , G7 , G8 , G9, GIO , Gil  CHA.L7  06 

,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23  CHA1705 

,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,C35  CHA1706 

REALX8  NG,MU2  CHA1707 

COMMON  /GRADS/ OUDXIL , DPDXIL , DGDXIL , DRDXI L , DZDXI L , DSDXIL ,  CHAI708 

L  DUDXIR, DPDXIR, DGDXIR, DRDXIR. DZDXIR, DSDXIR  CHA1709 

COMMON  /STEPO/UL,PL,ROL,GL,UR,PR,ROR.GR,USTAR,PSTAR,  CHA1710 

R5TARL,RSTARR,GSTARl,GSTARR,  CHA1711 

CL,CR,CSTARL,CSTARR, SL  ,  SR,  WL  ,  WR, UWC6  )  CHA1712 

, LAMDAL , LAMDAR, RATEL , RATERr TEMPI , TEMPR, TEMPSL , TEMPSR  CHA1713 
,  ZL ,  ZR, ZSTARL  ,  Z5TARR ,  NFLUX ,  HEL  EiiL  ,  HEL  EMR  CHA1714 

REALX8  LAMDAL, LAMDAR  CHAI715 

LOGICAL  HELEML , HELEMR  CHA1716 

COMMON  /STEP1/DUIDT, DPIDT, DGIDTL , DGIDTR, DRIDTL , DRIDTR  CHA1717 

2  , ASTARL , ASTARR, LAMDSL , LAMDSR, DSDAL , DSDAR , DZDAL , DZDAR  CHA1718 

3  , RAT , SH  CHA1715 

4  , BETACL , BETACR, DSDASL , DSDASR, DZDA3L , DZDASR  CHA1720 

REAL*8  LAMDSL, LAMDSR, DSDAL, DSDAR, DZDAL, DZDAR  CHA1721 

COMMON/ DETO/Q DET , PCJDET , RCJ  DET , UCJ  DET , DCJDET , PODET, ROODET ,  CHA17  22 

1  RATE, TFMPC  CHA1723 

COMMON  /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN,ZXN  CHAI724 

1  , GIH  CHA1725 

2  , FIH4 , ZMDOTL , ZMDQTR  CHA1726 

REALMS  LAMDAO  CHAI727 

Cxxx«xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxc.HA1728 

C  RO,U, P,Z  AND  THEIR  (XI, T)  DERIVATIVES  AT  EULFRIAN  POINT  X=X(I).  CHAI729 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHA1730 
DT2=DT/2 . DO  CHA173I 
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non  non  ooo 


CHAROEP 


FORTRAN  A1 


PAOE 


00  TO  C1,2,3,4,5,$),NFLUX 
1  CONTINUE 

NFLUX*1.  LINE  X«0  IS  TO  THE  LErr  OF  LEFT  HAVE. 

UX*UL 
PX*PL 
R0XuR0L 
ZX-Zt 
OX*OL 

DUDXIX'DUDXIi. 

DPDXIX«DPDXIL 
0R0XIX*DRDXIL 
DZDXIX-DZDXIL 
DUDTX--DPDXII 
DR0DTX»-R0L*«2KDUDXIL 
DPDTX— 0L**2*DUDXIL 
DRQDTX*DRODTX-KAT*ROL*UL 
DPDTX»DR0DTX*CLXK2 
DZDTX®0 . 

00  TO  9 
CONTINUE 

NFLUX»«.  LINE  X*0  IS  TO  THE  RIGHT  OF  RIGHT  HAVE. 

UX*UR 
PX«PR 
ROX-ROR 
ZX«ZR 
GX»GR 

DUDXIX*DUDXIR 
DPDXIX*DPDXIR 
DRDXIX»DRDXIR 
DZDXIX»DZBXIR 
DUDTX=*~DPDXIR 
DPDTX '~0RX*2*DUDXIR 
DRGDTX3~nORxx2*DUDXIR 
DRODTX=DRODTX-RATXROR*UR 
0PDTX=DRCDTX*CR*K2 
DZDTX=0. 

GO  TO  9 
CONTINUE 


N FI UX-2 .  SONIC  CASE  (LEFT). 

BETA0=(MU2X(UL/CL+G7 ) )*X(I .D0/MU2) 

SQBO=DSQRT ( BETAO ) 

A1=DUDXIL+DPDXIL/Gl 

A0=A1-(CL/(G15*SL) )*DSDXIL*(BETA0X*C5-1 .DO) 

EVER1  =  -C (GAMA- 1 .DO)XUL+2.DOXCL)X(BETAOXXG13-1.DO )/( GAMA-5 . DO ) 
EVER2=4 . D0XCLX(BETA0X*G14-1 . DO )/< 3 . D0*GAMA-5 . DO) 
EVER=(EVER1+EVER2)*RAT/R0l 
AO  =  (A(l+EVER) 

DPDAX”GLXBETA0*A0 
C0=MU2XCUL+G7*CL) 

IFtCO.LT.O.)  CALL  SOFC’FLUXE  2.  CO  NEGATIVE.’) 

UX=CO 

P.OX  =  GLXBETAO/CO 
ZX=ZL 

PX=RCIX*C0»X2/GAMA 
GX=R0XXC0 

DPDAX=DPDAX+RAT*UXXCOXSOBO 
DUDBX=-CLXBETA0XX(-1 . D0/G4)/G4 
DPDBX=PLXBETA0*XMU2/G6 
DRODBX=ROL*BETAOX*(-MU2)/G4 
DSDAX=SQBO*DSDAL 
DZDAX=SQBO*DZDAL 

DRQDAX=DPDAX/COXX2-(ROX/(GAMAX5L))XDSDAX 
DUDAX=AO 

DGDAX=G . 5DO*GAMA*(PX*DRODAX+ROX*DPDAX)/GX 
GO  TO  9 
5  CONTINUE 


CHA1732 

CMA1733 

CHAI734 

CHA1735 

CHA1736 

CHA1737 

CHA173R 

CKA1739 

CHA1740 

CHAI741 

CHA1742 

CHA1743 

CHA1744 

CHAI745 

CHA1746 

CHA1747 

CHA1748 

CHA17**9 

CHA1750 

CHA1751 

CHAI752 

CHA1753 

CKA1754 

CHA1755 

CHA1756 

CHA1757 

CHA1758 

CHA1759 

CHA1V6Q 

CHA1761 

CHA1762 

CHA1763 

CHA17&4 

CHA1765 

CHA1766 

CHA1767 

CMA17o8 

CHAI/69 

CHA1770 

CHA1771 

CHAI772 

CHA1773 

CHA1774 

CHA1775 

CHA1776 

CHA1777 

CHA1778 

CHA1779 

CHA17SG 

CHAI781 

CHA1782 

CHA1783 

CHAI784 

CHAI785 

CHA1786 

CHAI787 

CHA178S 

CHA1789 

CHA1790 

CHA1791 

CHAI792 

CHAI793 

CHA1794 

CHA1793 

CHA1796 

CHAI797 

CHAI798 

CHA1799 

CHA1800 

CHA1801 

CHA1S02 

CHA1893 
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non  non 


CNAROEP  FORTRAN  M 


NFLUX-5.  SONIC  CASE  (RIGHT). 

BET A0»  <MU2*( -UR/CR+G7 ) > «*( 1 . DC/MU2 ) 

S9B0*DSQRT(BETA0) 

A1*DUDXIR-DPDXIR/9R 

AO»Al+CCR/CG15*SRJ)KDSDXIR*(BETAOX*G5-l.DO) 
EVER1*(-(OAHA-1.DO)XUR+2.DOXCR)XCBCTAOXXG13-I.DO)/(OAMA-3.DO) 
EVER2*-4 . DOXCRXC  DETA0XXG14-1 . DO )/<  3 . DOKQAMA-3 . DO) 
EVER*(EVERl+EVER2)*RAT/R0R 
AO-CAO+EVER) 

DPDAX»~GR*BETA0*A0 
C0*MU2*( -UR+G7XCR) 

IF(C0  . LT.O. )  CALL  SOFCFLUXE  5.  CD  NEGATIVE.') 

UX»-C0 

ROX>GRMBETAO/CO 
ZX*ZR 

PX»ROX*CO**2/GAMA 
GX«ROX*CO 

DPOAX-DPDAX-RATXUXXCO»DSQRT<BETAO) 

DUDBX*CRXBETAflXX(  -1 . D0/G4 )/G4 
DPD8X*PRXBETA0K*MU2/G6 
DRODBX«ROR*BETAO**(-MU2)/G4 
DSDAX»SQBO*DSDAR 
DZDAX«3QB0*DZDAR 

DR0DAX»DPDAX/C3**2-(R0X/<GAMA*SR)>*DSDAX 
DUDAX*AO 

DGDAX°0 . 5DOXUAMAX(PXXDRODAX+ROXXDPDAX)/GX 
GO  TO  9 
CONTINUE 

NFLUX«3.  LINE  X*'Q  IS  BETWEEN  THE  LEFT  WAVE  AND  THE  CONTACT. 

UX*USTAR 
PX=PSTAR 
ROX=RSTARL 
ZX=ZL 
GX=GSTARL 

DUDXIX=-DPIDT/GSTARL*X2 
DPDXIX*-DUIDT 

DUDXIX*DUDXIX-RATXU5TAR/RSTARL 
DZDXIX=DZDXIL 
DZDTX=0. 

IF  ( . NOT . HELEML )  GO  TO  32 

31  CONTINUE 
C  LEFT  SHOCK. 

DRDXIX=(RSTARL/WL)X*2*(3. DOXDUIDT 

1  +DPIDTX(1 .D0+3.D0*CWL/GSTARL)XX2)/NL 

2  +DUDXIL*WL*C(GL/WL)**2+3.D0)+3.D0*DPDXIL 

3  +DRDXIL*CWL/R0L)*X2) 

EVER1=UL*RST/\RLXX2xRAT*<  (GL/WL)**2+1 . DO )/( ROLXWL ) 

EVER2  =  2 .  DO*R'SrARL*USTAR*RAT/WL 

DRDXIX=DRDXIX+EVER1+EVER2 

DR0DTX=-DUDXIX*R0X**2 

GO  TO  33 

32  CONTINUE 
BETA--G3TARL/GL 
SQB=DSQRT(BETA) 

DPDA=ASTARLXGSTARL 

DPDA=GSTARL*(ASTARL+RATXUSTARXCSTARL/CGL*  SQB) ) 
G41=l.D0/G4+0.5Du 

DRODA  =  (DRDXIl-DPDXIL/<CL*CD)  XBETAX*G41  +  DPDA/'CC3TARLX*2) 

DRDXIX=  DR0DA/SQB+DPIDT/(GSTARt*CST4RL**2) 

DR00A  =  DPDA/CSTARLXX2-(RSTARL/CGAMAXSU)XDSDASL 
DROD  TXS-DUDXIXXR0XX*2 
DRDXIX=DRODA/SQB+DRODTX/GSTARL 

33  CONTINUE 
DUDTX=DUI DT 
DPDTX=DPIDT 
GO  TO  9 

4  CONTINUE 
C 


CHA1304 

CHA1805 

CHA1806 

CHA1307 

CHA1808 

CHA1809 

CHA1810 

CHA1811 

CHA1812 

CHA1813 

CHA1814 

CHA1815 

CHA1816 

CHA1817 

CHA1818 

CHA1819 

CHA1320 

CHA182I 

CHA1822 

CHA1823 

CHA1824 

CHA1825 

CHA1826 

CHA1827 

CHA1628 

CHA1K29 

CHA1830 

CHA1831 

CHA1332 

CHA1833 

CH41834 

CHA1835 

CHA1836 

CHA1837 

CHA1838 

CHA1839 

CHA1840 

CHA1341 

CHA1842 

CHA1843 

CHAI844 

CHA1845 

CHAj.846 

CNA1347 

THA1848 

CHA1849 

CHA1850 

CHA1851 

CHA.U52 

CHAI853 

CHA1854 

CHA1S55 

CHAI856 

CHA1857 

CHA1858 

CHA] 859 

CHA1360 

CHA1861 

CHA1862 

CHA1863 

CHA1864 

CHA1865 

CHA1S66 

CHA1867 

CHA186S 

CHA1869 

CHA187U 

CMA1871 

CHAI372 

CHAI  73 

CHA1374 

CHA1875 
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CHARGE? 


FORTRAN  A1 


C  NFIUX"4,  LIRE  <-0  IS  BETWEEN  THE  CONTACT  AND  THE  RIGHT  WAVE. 
C 

DPDXIX—DUIDT 

UX’USTAK 

PX»rSTAR 

ROX“RSTARR 

2X«ZR 

GX*GSTARR 

0UDXIX«-”0PIDT/GSTARR*»2 

DUDXIXUDUDXIX~RATRUSTAR/RSTARR 

DPDXIX*-DUIDT 

DZDXIX“DZDXIL 

DZDTX-0. 

IF  ( . NOT.HELEMR)  00  TO  42 

41  CONTINUE 
C  RIGHT  SHOCK 

DRDXIX3  <  R5TARR/WiT  )  **2*U  3 ,  KDU I OT 

1  -DPIDTHC 1 . DQ+S . DUKU4R/QSTARR ) x«2 )/WR 

2  -DUDXTR*HRK(  CGR/WR)x*2+3 . DO )+3 . DQXDPDXIR 

3  +DRDXIR*(WR/RQR)**2> 

EVER1 ‘iUR*RSTARRK*2*RAT*{ ( GR/WR) **2+1 . DO  V ( RGR*NR) 

EVER2*2  DOXRSTARRKUSTARXRAT/HR 
DRDXXX-DRDXIX-  EVER1-EVER2 
DRODTX  -DUDXIXKR0XXK2 
GO  TO  43 

42  CONTINUE 

C  RIGHT  RAREFACTION 
BETA-G3TARR/0R 
SQB*DS«RT(3ETA) 

DPPA*-ASTARR*G5T  arr 

DPDA«-OSTARR*( ASTARR+RAT*USTAR*CSTARk/(GR*  SOB)  ) 

041*1. DO/G4+0.5DO 

DRQDA*(DRDXIR-DPDXIR/(CR*CR) )  *BETA**G41+DPDA/'CSTARR**Z) 

DRDXIX*  DR0DA/SQB-DPJDT/(G$TARR*CSTARR*K2) 

DR0B«*DPDA/CSTARR**2-(RSTARR/<GAMA*SR))XD3DASR 

DR0DTX=-DUDXIX*R0X**2 

DRDXIX=DRODA/SqB-DRGDTX/GSTARP 

43  CONTINUE 
DUDTX-DUinT 
DPDTX=DPIDT 
30  TO  9 

9  CONTINUE 


CHA1876 

CHAX877 

CHA1878 

CHA1879 

CHA1880 

CHA1A81 

0HA1882 

CKAJ883 

CHA1884 

CHA1885 

CHA1886 

CHA1887 

CHA1888 

CHA1889 

CHAI890 

CHA1891 

CHA1892 

CHAI893 

CHA1894 

CHA1895 

CHA1896 

CHA1897 

CHA1S98 

CHA1899 

CHA1900 

CHA1901 

CHA1902 

CHA1903 

CHA1904 

CHA1905 

CHA1906 

CHA1907 

CHA1903 

CHA1909 

CHA1910 

CHA1911 

CHA1912 

CHA1913 

CHA1914 

CHA1915 

CHA1916 

CHA1917 


C*XXXXX**#*X*XKXX#***X**XX*:fXXXX*X***3<X*X*XXX*XXXXXX*XXXXXXXXK**KXXXXXXXCHA19:i8 


C  FLUXES  CENTERED  AT  TIME  T<N+l/2)  AT  EULERIAN  POINT  X=X(I)  CHA19I9 

CXXXXXKXX.'tXXXXXXXXXXXXXXXX*<XXXXX:'<XXXXXXX:XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCHAl?20 


FI1=RQX*UX 

FI2=RUXXUXXJ'.2+?X 

FI2-»FJ:2-PX 

FI3=UXXIG12XPX^0.5D0XR0XXUXXX2) 

FI4=ZXxR0XxUX 

FI3aCI3+QDETxFI4 

ROUOO=ROX*UX 

GO  T0C1C, 20, 30,40,50,60),  NFLUX 
10  CONTINUE 
60  CONTINUE 

DFaXIl=DRDXIX*UX+ROXXDUDXIX 
DFDXI2=DRDXIXxUXXX2+2 .  DO^RQXXUXXDUDXIX+DFDXJX 
DFDXI2-DFDXI2-UFDXIX 

DFDXI5=DUDXIXXCG12xPX+0 .5D0XROXXUXXX2) 

1  +UXX(G12xflPDXIX+C  . 5D0XDRDXIXxUXxx2+R0XXUXXDUDXIX) 

DFDXI4=ZXXDFDXII+R0XXUXXDZDXIX 
DFDXI3  =  DFDXI  3+QI)ETxDf  DXI4 
DFIDTl=DRODTX#UX+RaX*DUDTX 
DFIDT2=DRQDTXXUXXX2+2.D0*RQXXUXXDUDTX+DPDTX 
DFIDT2=DFIDT2~DPDTX 

DFIDT3=DUDTXX(G12*PX+0 . 5D0*R0X*UX*#2 ) 

1  +UX*(G12XDPUTX:  Q.5DO*DRODT/XUX<*2+ROX*UX*DUD7X) 

DFIDT4-7.XXDFinTl+R0X*ZX*DZDTX 
DFIDT3=DFI DT3+QDE7*DFIDT  4 
FI!)OT1=-ROUOOXPFDXI1  +  DFIDTI 
FIPOT2=-ROUOOXDFDXI2+DFIDT2 
FIDOT3--ROUOQXDFDXI3+DFIDTS 


CHA1921 

CHAI922 

CHA1923 

CHA1924 

CHA1925 

CHA1926 

CHA1927 

CHA1928 

CHA1929 

CHA1930 

CHA1931 

CHA1932 

CHAI933 

CHA1934 

CHA1935 

CHA1936 

CHA1937 

CKA1938 

CHA1939 

CHA1940 

CHA194I 

CHA1942 

CHA1943 

CHA1944 

CHA1945 

CHA1946 

CHA1947 
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me.  CHAROEP  FORTRAN  A1 


FID0T4—R0U0Q*BFDXI1+DFIi9T4 

UX30T*-ROUOOKRUDX7.X+DUDTX 

PXDOT«*-ROUOOKOPmX+DPDTX 

ROXD0T«~RQUOO*DRDXJX+DRODTX 

ZXOOT*-RmO*DZDXIX+DZDTX 

FIN3«FU+DT2*FID0T1 

FIM2“FI2<DT2*FID0T2 

GIH»PX+3T2*PXD0T 

FIH3°FI3AErV2«FIt'0T3 

FIH4*FI4«-DT2XFID014 

UXN«UX+DT*UXDGT 

PXN-PX+DTXPXDOT 

ROXN-SOX+DTXROXDOT 

2XN-2X+DTNZXD0T 

IFCZXN.IT.O.)  ZXNa0 . 

(JO  TO  90 
20  CONTINUE 

EVO«0UB5«RT(BETAO) 

201  CONTINUE 

DFIDA1»DR0DAX*UX+R0X*DUDAX 
DFIDA2*DR0DAX*UXRX2+2 . D0*ROXXUX*DUDAX+DPDAX 
DFIOA2*OFXDA2-DPMX 
DFIDA3*DUDAX*(0J,2XPX+Q.5DO*ROXXUXH*2) 

1  +UXR(O12*DPDAX+0 . 5nC*0R0DAX*UXXX2+R0X*UXX&*JDAX) 

DFIDA4*ZX«DFIDA1+R0X*UXXDZDAX 
F3DOT1«~EVO*DFJDA1 
FID0T2*"-EV0X’)FIDA2 
FIDGT3*~EV0*DFIDA3 
F3DOT4—EVOXDFIOA4 
FTH1*FI J  +BT2XFIG3T1 
FIH2*FI2+DT2*FIDUT2 
FIH3«FI3+BT?*FXDOT3 
FIH4*FI4+DT2*FI!inT4 
GA-3GDAX 

IF(NFl.UX.EQ..5)OA»-GA 

DROUA«UX*DRODAX+ROXXDUDAX 

BETAPR*0.5D0xDSQRT(BETA0)«CCfA-UROUA) 

FIH2  =  FIK2--DPDBXJl3ETAPRXDr2 
JXDOT =-EV0XDU0AX+BETAPRXBUDSX 
rXDOT-’-EVOXDPDAX+EETAPRxDPDBX 
GIH=PX+DT2*PXDQT 

ROXDOT3-£Va*CRODAX+6ETAPR*DRODBX 

ZXDOT=-EVO*DZDAX 

UXN3UX+DT*UXDOT 

PXN=PX+DTxPXDOT 

R0XN3R0X+DT*R0XEQ1 

ZXN=ZX+DT*ROXDOT 

IFCZXN.LT.O. )  ZXN-0. 

GO  TO  90 
50  CONTINUE 

EVQ=“GRXDSQRT ( BETAO ) 

GO  TO  201 
50  CONTINUE 
40  CONTINUE 
GO  TO  60 
90  CONTINUE 
RETURN 
END 
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CHA1951 
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FAKE  FACTION  BRANCH 
SHOCK  £>P*/4CH 


Figure  A-2.  Intersection  of  Right  and  Left  Adiabats  for  Solving  Riemann  Problem 


54 


WIVI1 


APPENDIX  B.  Code  for  Re- Normalizing  the  Air  Impulse 


3 


1  IMPLICIT  REAL*8(A-H,0-Z)  REN00010 

C  CODE  RENORM  —  C  TRANSFORMATION  OF  TOTAL  REFLECTED  IMPULSE  FROM  REN00020 
C  BAKER'S  CHART  TO  SPACE-NORMALIZED  VALUES.  REN00030 

C  DATA  FROM  FIO.  6.3  (SUPPLEMENT)  IN  BAKER'S  BOOK  "EXPLOSIONS  IN  AIR"  REN00040 
Z  REALK4  RB, IB, RS, IS, ISBARE  REN00050 

3  DIMENSION  RB<21 ) , IB<21 )  REN00060 

4  DIMENSION  RS(21 ) , IS(21 ) , ISBAREC21 )  REN00070' 

5  DATA  RB/.05, .06, .07/ .08, .09, .1, .2, .3, .4, .5, .6, .7, .8, .9,1.,  REN00080 

1  2. ,3. >4. ,5. ,6. ,7./  REN00090 

6  DATA  IB/4.4,3.06,2.30,1.83,1.50,1.27,.457,.293,.221,.178,.149,  REN00100 

1  .128, .113, .099, .0885, .0376, .0236, .0173, .0136, .0113, .0095/  REN00110* 

7  PAI«4.D0*DATAN(1.D0)  REN00120 

8  G*1 . 4D0  REN00130 

9  PA*0.ID0  REN00140 

10  RHOA«1.3DO  REN00150 

11  RHOO-18QO.DU  REN00160 

12  90*4 . DO  REN00170 

13  BETA»DSQRT(RHOA/RHOQ)*<PA/<RHOO*QO))**<1.D5/&.DO)  REN00180 

14  OOREM»(3.DO/DS9RT(2.DOXG))*C4.DOXPAI/S.DO)X#a.DO/3.DO)  REN00190 

15  BETA-BETAXOOREM  REN00200 

16  DELTA*(  (4.DOXPAI/3.DO)H(RHOOXQO/PA)  )XX( 1 . DO/3 . DO)  REN00210 

17  PRINT  11,  BETA, DELTA  REN00220 

18  11  FURMAT(/1X, 'RESULTS  WITH  BETA, DELTA* ' , 2D16 .  7'/  REN00230 

1  IX,'  N','  RB  ','  IB  • , 2X,  REN00240 

2  '  RS  ','  IS  • , 2X,  '  ISBARE  V)  REN00250 

REN00260 

19  DO  1  Nal , 21  REN00270 

20  RS(N)*RB<N)XDELTA  REN00280 

21  IS(N)-IB(N)XBETA  REN00290 

22  ISBARE(N)*1.D0/RS(N)XX2  REN00300 

23  PRINT  2,  N,RB(N),IB(N),R5(N),IS(N),ISBARE(N)  REN00310 

24  2  FORMAT (IX, 14, 2E12 . 4 , 2X,2E12 . 4, 2X, E12 . 4)  REN00320 

25  1  CONTINUE  REN00330 

26  END  REN00340 


RESULTS  WITH  BETA 

l,  DELTA*  0, 

.12041630-01 

0.6706157D+02 

N 

RB 

IB 

P.S 

IS 

ISBARE 

I 

0.5000E-01 

0 . 440UE+0.1 

0.3353E+Q1 

0 . 5298E-01 

0.8894E-01 

2 

0.6000E-01 

0.3060E+01 

Q  .4024E+01 

0.3685E-01 

0 . 6 177E-01 

3 

0 .7000E-01 

0 .2300E+01 

0.4694E+01 

0 .2770E-01 

0 .4538E-01 

4 

0 . 8Q00E-Q1 

0 . 133GE+01 

0 . 5365E+0i 

0.2204E-01 

0 .3474E-01 

5 

0.9000E-31 

0.1500E+01 

0 .6036  E+01 

0. 1806E-01 

0 . 27  45E-01 

6 

0 .lOOOE+OU 

0. 1270E+G1 

0 . 67 06E+01 

0 . 1529E-01 

0 .2224E-01 

7 

0.2000E+00 

0 . 4570E+00 

C  -1341E+02 

0 . 5503E-02 

0.5559E-02 

8 

0 .3030E+0Q 

0 .2930E+00 

0 .2012E+02 

0.3528E-02 

0.2471E-02 

9 

0 . 4000E+30 

0.2210E+Q0 

0.26S2E+02 

0.2661E-02 

0 . 1390E-02 

10 

0.5000E+00 

0.178QE+QU 

0 . 3353E+02 

0.2143E-02 

0.8894E-03 

11 

0.6000E+00 

0 . 1490E+00 

0. 4024E+02 

0 . 1794E-02 

0.6177E-03 

12 

0.7000E+00 

0.1280E<-G0 

0 .  **6  94E+C2 
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